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CONSTITUTIVE  MODELING  OF  ROCKS 
WITH  INTERNAL  CRACKS  AND  PORES 


1.0  INTRODUCTION 


The  brittle  deformation  processes  in  materials  such  as  rocks,  concretes,  ceramics,  etc.,  have 
lately  attracted  a  great  deal  of  attention.  On  the  microscale  this  class  of  deformation  processes 
is  characterized  by  the  cooperative  evolution  of  a  large  number  of  crack-like  microdefects  of 
irregular  geometry.  The  complexity  of  the  problem  is  further  exacerbated  by  the  dependence 
of  the  crack  growth  and  the  stifihess  on  the  sign  of  normal  stresses.  Early  constitutive  models 
were  based  on  modifications  of  the  venerable  plasticity  theory  originally  intended  to  predict  the 
behavior  of  ductile  metals.  The  inability  of  these  models  to  predict  the  brittle  response  of  rocks 
stimulated  efforts  to  examine  the  underlying  micromechanical  processes.  Constitutive  analy¬ 
sis  which  captures  salient  aspects  of  the  xmcromechanical  processes  characteristic  of  particulu 
materials  and  applications,  which  is  still  reasonably  simple,  is,  therefore,  highly  desirable. 

The  early  micromechanical  models  addressed  primarily  the  dependence  of  elastic  moduli  on 
the  microdefect  density.  Budiansky  and  O’CoxmeU  (1976)  considered  randomly  distributed 
flat  cracks  and  applied  the  self-consistent  method  to  obtain  the  effective  elastic  properties  of 
the  overall  isotropic  response.  Li  this  analysis  all  cracks  were  supposed  to  be  open  during 
loading.  In  the  case  when  some  cracks  dose  or  undergo  frictional  sliding,  the  overall  response 
becomes  anisotropic  and  load-path  dependent  (Horii  and  Nemu.t-Nasser,  1983).  Sammis  and 
Ashby  (1986)  devdoped  an  approximate  theory  to  predict  the  failnre  of  brittle  porous  solids, 
loaded  in  compression,  during  which  cracks  grow  from  the  surfaces  of  the  preexisting  pores  or 
vacandes.  Nemat-Nasser  and  Obata  (1988),  in  thdr  micromechanical  modeling  of  the  inelastic 
response  of  brittle  materials  such  as  compact  rocks,  concrete  and  some  ceramics,  used  the  sliding 
crack  mechanism  as  the  dominant  source  of  indastidty.  According  to  this  modd  the  frictional 
diding  of  preexisting  cracks  leads  to  the  formation  (ff  tension  wing-cracks  (see  also  Kachanov, 
1982,  and  Ashby  and  Hallam,  1986).  A  dilute  distribution  of  preexisting  cracks  was  assumed, 
neglecting  interaction  among  neighboring  flaws. 

Other  mechanisms  of  microcracking  in  rocks,  such  as  dastic  mismatch  and  bending  mechanisms, 
were  also  suggested.  For  heterogeneous  materials,  the  local  tensile  stresses  appear  at  the  interface 
of  two  elastically  mismatched  materials.  Tensile  stresses  resulting  from  the  unequal  lateral 
expansion  of  the  mismatched  materials  are  often  suffident  to  cause  nudeation  and  subsequent 
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propagation  of  a  microcrack  (Kemeny  and  Cook,  1991).  According  to  the  bending  model, 
the  tensile  stress  needed  to  trigger  microcrack  growth  occurs  as  a  result  of  bending  a  soft 
and  elongated  particle  spanning  two  harder  inclusions.  Extensive  summaries  and  reviews  on 
microcracks  in  rocks  are  given  by  Kranz  (1983),  Zheng  (1989),  and  others. 

The  behavior  of  brittle  rocks  with  inferior  tensile  strength  depends  on  the  mode  and  stability  of 
the  crack  growth.  Thus,  the  response  of  a  brittle  rock  subjected  to  compressive  loading  strongly 
depends  on  the  lateral  confinement.  An  unconiined  specimen  containing  a  large  number  of  flaws 
fails  by  axial  splitting  or  slabbing  (Horii  and  Nemat-Nasser,  1986).  Final  failure  occurs  as  a 
result  of  the  unstable  growth  of  a  single  crack  of  preferential  geometry  at  a  relatively  small 
microcrack  density.  Thus,  the  microcrack  interaction  is  relegated  to  a  second-order  effect.  As 
the  confinement  is  increased,  axial  splitting  is  suppressed  and  a  typical  rock  specimen  fails  by 
formation  of  a  narrow  re^on  of  high  crack  density  (fault  or  "shear  band”).  Finally,  at  large  levels 
of  lateral  confinement,  homogeneous  microcracking  prevails  throughout  the  sample,  resulting  in 
a  quasi-ductile  overall  response.  The  microcrack  density  (especially  within  the  "shear  band” 
emerging  near  the  apex  of  the  force-displacement  curve)  is  in  this  case  much  more  substantial, 
rendering  microcrack  interaction  not  only  important  but  even  the  dominant  factor  of  failure.  The 
degree  of  confinement  at  which  the  mode  of  failure  changes  is  referred  to  in  rock  mechsmics  as 
the  brittle-to-ductile  transition.  This  phenomenon  was  studied  by  many  authors,  most  recently 
by  Hirth  and  Tnllis  (1989),  Wong  (1990),  and  Zhang,  et  aL  (1990a,b). 

The  micromechanical  studies  reported  in  the  literature  have  demoxutrated  all  of  the  advantages 
and  disadvantages  of  micromechanical  models.  The  relative  absence  cf  ambiguity,  vlirect  identi¬ 
fication  of  material  parameters,  and  conceptual  clarity  of  micromechanical  models  are  accompa¬ 
nied  by  a  lack  of  computational  efficiency  and  limitations  to  simple  geometries  and  homogenous 
states  of  stress.  By  its  very  nature  the  imcromechanical  determination  of  the  state  of  damage  at 
a  material  pdnt  of  an  effective  continuum  involves  compilation  of  records  defining  the  growth  of 
each  microcrack  within  the  corresponding  representative  volume  element  of  the  actual  material. 
The  attendant  bookkeeping  and  averaging  represents  a  formidable  problem  which  may  be  re¬ 
dundant  since  it  is  not  immediatdy  known  how  detailed  the  information  sought  must  be.  Thus, 
as  stated  by  Bice  (1975),  a  direct  micromechanical  prediction  of  material  response  is  unlikdy 
to  displace  the  application  of  phenomenolofpcal  and  less  rigorously  based  structure-parameter 
modds. 

Unfortnnatdy,  the  state  of  the  art  of  phenomenolopcal  modeling  is  far  from  satisfactory.  In  the 
wake  of  the  original  Kachanov  modd  ( 1958),  in  which  damage  was  measui-ed  by  a  scalar,  different 
and  contradictory  choices  of  other  damage  variables  were  promoted  as  proper  characterizations 
of  the  deterioration  of  the  material  properties.  Vectorial,  second-order,  and  fourth-order  tensor 
representations  of  damage  were  proposed  to  capture  the  damage-induced  anisotropy  (Vakulenko 
and  Kachanov,  1971,  Dragon  and  Mroz,  1979,  Kachanov,  1980,  Krajcinovic  and  Fonseka,  1981, 
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Talreja,  1985,  Cottin,  1985,  Murakami,  1988,  Ju,  1989,  Kachanov,  1992,  etc.).  The  phenomeno¬ 
logical  features  of  damage  mechanics  theory,  such  as  damage  surface  and  damage  rule,  analogous 
to  yield  surface  and  flow  rule  of  the  more  developed  phenomenological  metal  plasticity  theory, 
are  still  in  the  early  stages  of  development  and  are  the  subject  of  ongoing  active  research  (Simo 
and  Ju,  1987,  Chow  and  Wang,  1988,  Krajcinovic,  1989,  Chaboche,  1992,  etc.).  Morover,  the 
relation  of  these  models  to  the  physics  of  the  phenomenon  is  often  less  than  obvious. 

A  third  class  of  models  makes  an  attempt  to  combine  the  desirable  attributes  of  micromechan¬ 
ical  and  phenomenological  theories.  Rudnicki  and  Rice  (1975)  treat  the  onset  of  rupture  in 
brittle  rock  masses  as  a  constitutive  instability.  Their  macroscopic  constitute equations  are 
the  pressure-dependent  generalization  of  the  Prandtl-Reuss  equations  of  metal  plasticity.  In¬ 
elastic  response,  dilatant  in  nature,  was  considered  to  be  predominately  a  consequence  of  the 
frictional  sliding  on  microcrack  surfaces,  their  uplifting  over  asperities,  and  tensile  cracking  into 
the  crack  wings  (kinks).  Nemat-Nasser  and  Shokooh  (1980)  further  developed  these  constitu¬ 
tive  equations  to  account  for  inelastic  volume  changes  and  pressure  sensitivity,  important  in 
geotechnical  materials  such  as  cohesionless  sands  and  cohesive  soils.  Ortiz  (1985)  considered 
a  constitutive  model  for  the  inelastic  behavior  of  concrete,  by  viewing  it  as  a  mixture  of  two 
phases:  mortar  and  aggregate.  Simple  modds  were  used  to  describe  the  individual  responses 
of  mortar  and  aggregate,  and  mixture  theory  was  implemented  to  obtain  the  overall  composite 
response.  None  of  the  existing  models,  however,  satisfies  the  contradicting  requirements  of  rigor 
and  simplicity.  Moreover,  none  of  these  models  was  used  successfully,  if  at  all,  for  large-scale 
computations  or  for  a  general  case  of  loading. 

The  objective  of  the  present  study  is  neither  to  recapitulate  the  current  state  of  micromechan¬ 
ical  and  phenomenological  modeling  nor  to  reconcile  all  or  even  most  existing  theories  singling 
out  those  that  do  not  satisfy  all  conditions.  Instead,  the  ultimate  objective  is  to  formulate  a 
microscopically-inspired  continuum  model  which  will  be  suitable  for  large-scale  computatioiu. 
While  retaining  as  much  generality  as  possible,  the  proposed  model  will  specifically  address  the 
brittle  defer  - a  processes  in  porous  sandstone  and  limestone  rocks.  Ductile  deformation  will 
be  conndexeu  ir  ,  following  study.  In  view  of  the  many  competing  models,  it  seons  reasonable  to 
revisit  the  problem  from  its  very  start  and  examine  the  approximations  related  to  the  representa¬ 
tion  of  the  existing  distribution  of  microcracks  by  an  appropriate  tensor  measure.  Additionally, 
the  present  study  will  reexamine  the  mathematical  structure  of  the  effective  compliances  and 
sti&esses  demonstrating  that  rigor  can  be  preserved  without  attendant  sacrifices  in  simplid- 
ty.  At  this  pdnt  a  second-order  tensor  representation  of  damage  will  be  considered  and  the 
damage-induced  anisotropy  will  be  approximated  by  orthotropy.  A  more  complete  discussion  of 
the  problem,  within  the  framework  formulated  herein,  including  considerations  o' other  types  of 
anisotropy  and  extensions  to  rate  theories,  will  be  presented  subsequently. 
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2.0  MICROMECHANICS  OF  BRITTLE  DEFORMATION 
PROCESSES  IN  POROUS  ROCKS 


The  word  rock  refers  to  a  large  class  of  naturally  formed  solid  materials  of  different  structure  and 
chemical  composition.  Near  the  exposed  surface,  most  rocks  are  weathered  and  fragmented.  The 
stresses  to  which  the  rocks  are  exposed  in  engineering  application  are  often  inferior  to  stresses 
developed  during  their  formation  or  tectonic  movements.  As  a  result,  a  typical  rock  contains 
crack-like  defects  on  a  range  of  scales.  From  the  geolo^cal  and  seismic  viewpoint  the  most 
important  defects  are  joints  and  faults  with  lengths  measured  in  tens  and  hundreds  of  miles. 
Our  current  interest,  however,  is  centered  on  studying  the  influence  of  relatively  small  defects  on 
the  mechanical  response  of  rock  specimens  which  do  not  contain  joints,  tested  in  the  laboratory. 

The  mechanical  response  and  strength  of  a  rock  is  to  a  large  degree  determined  by  its  microstruc¬ 
ture  (or  fabric).  The  fabric  of  a  typical  rock  encompasses  crystal  aggregates  joined  together  by 
some  cementitious  material.  Preferred  crystal  orientations,  bedding  planes,  foliation,  schistosity, 
microfissures,  and  porosity,  common  to  rocks,  are  the  principal  reasons  for  preferred  directions, 
anisotropy,  and  reduction  of  the  mechanical  strength  (see,  for  example,  Jaeger  and  Cook,  1976). 

As  a  consequence  of  their  inferior  tensile  strength  rocks  are  subjected  primarily  to  compressive 
loadings  in  en|^eering  applications.  The  average  (macro)  normal  stresses  are  seldom  if  ever 
tensile.  However,  the  local  stress  fluctuations,  attributable  to  the  inhomogeneities  in  the  rock 
fabric,  will,  in  general,  have  tensile  components  as  well.  These  local  tensile  stresses  are  in  most 
cases  considered  to  be  the  principal  reason  for  microcracking.  The  distribution  and  magnitude 
of  these  local  stresses  d^ends  on  the  shape  cd  the  inhomogenrity,  dastic  mismatch,  defect 
density,  etc.  Thus,  the  determination  of  the  local  stress  from  the  average  stress  requires  a 
micromechanical  study  based  on  the  actual  rock  fabric. 

Some  of  the  definitions  used  in  rock  mechanics  are  a  direct  consequence  of  the  idiosyncrasies  of 
the  deformation  process  in  rocks.  The  conventional  crystal  plasticity  mechanisms  are  operative 
in  low-porosity  crystalline  rocks  only  at  high  temperatures  (Evans,  et  oL^  1990).  In  all  other 
cases  the  inelastic  deformation  is  attributed  to  the  nudeation,  propagation,  and  coalescence  of 
microcracks  into  larger  clusters.  For  example,  in  porous  rocks  macroscopic  ductility  reflects 
distributed  grain-scale  crushing  and  nucrocracking  (Wong,  1990).  The  mode  and  stability  of 
microcracking  depends  on  the  ratio  of  the  hydrostatic  pressure  and  deviatoric  stress.  Thus,  the 
ductile  deformation  of  porous  rocks  is  dependent  on  the  first  invariant  of  the  macrostress  tensor 
(see  Byerlee  and  Brace,  1969,  Jones,  1980,  etc.).  The  dnctOity  of  the  deformation  process  is 
customarily  defined  by  the  character  of  the  functional  dq>endence  of  the  axial  strain  on  the 
hydrostatic  and  deviatoric  stresses. 

In  stress-controlled  tests  at  low  confinement  the  mode  of  failure  is  brittle.  The  onset  of  failure 
is  sudden  and  is  not  preceeded  by  significant  inelastic  strains.  In  contrast,  the  much  more 
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important  behavior  of  rock  specimens  in  confined  conditions  is  rich  in  det^s  and  possibilities.  In 
crustal  conditions  the  deformation  of  rocks  in  the  majority  of  cases  is  controlled  by  displacement 
constraints  (confinement). 

Deformation  is  referred  to  as  being  ductile  on  the  macroscale  when  the  increase  of  the  axial  strain 
in  a  specimen  subjected  to  constant  hydrostatic  stress  requires  increase  of  the  deviatoric  stress. 
On  the  microlevel  this  typically  means  that  the  microcrack  growth  is  stable,  i.e.,  that  the  crack 
can  grow  only  if  the  deviatoric  stress  is  increased.  In  contrast,  if  in  a  displacement-controlled 
test  the  deviatoric  stress  reaches  a  maximum  and  starts  declining  at  increasing  axial  strain  the 
response  is  referred  to  as  being  brittle.  The  corresponding  segment  of  the  force-displacement 
curve  is  labeled  as  “softening”.  The  softening  phenomenon  is  readily  explained  on  the  microscale 
by  the  unstable  growth  of  microcrack  lengths.  Depending  on  the  ratio  of  the  differential  stress 
to  the  effective  mean  stress  the  ultimate  failure  at  low  temperatures  can  occur  as  a  result  of  the 
localization  of  microcracks  into  a  "crack  (or  shear)  band”  or  by  cataclastic  flow  (Paterson,  1978, 
Horii  and  Nemat-Nasser,  1986).  The  transition  between  the  failure  emphasizing  inhomogeneous 
deformation  (localization)  and  the  failure  characterized  by  homogeneous  deformation  (cataclastic 
flow)  is  referred  to  as  the  brittle-to-ductile  transition  (Rutter  and  Hadizadeh,  1991,  Wong,  et 
al.,  1992,  etc.).  At  high  confinement  levels  the  deformation  is  homogeneous  and  the  ultimate 
(or  macro)  failure  takes  place  when  the  density  of  microdefects  reaches  critical  concentration. 

2.1.  EXPERIMENTAL  OBSERVATIONS 

Even  though  the  relation  between  the  mode  of  deformation  and  microcrack  evolution  has  been 
accepted  as  the  dominant  aspect  of  the  phenomenon,  the  experimental  data  available  in  the 
literature  are  primarily  related  to  macroscopic  observation  (Hustruhd  and  Robinson,  1973,  Pa¬ 
terson,  1978,  Felice,  et  oh,  1991).  The  present  ^cusnon  (based  on  the  Hterature)  is  confined  to 
the  low-temperature  compression  of  sandstone  rock  specimens  below  the  brittle-ductile  transi¬ 
tion.  In  a  series  of  publicatioiu,  Zhang,  et  aL  (1989, 1990a,b),  and  Wong  (1990)  summarized  an 
ambitions  mcperimental  program  considering  several  different  sandstone  rocks  using  both  nonde¬ 
structive  testing  methods  (Acoustic  Emisuons  •  AE)  and  Scanning  Electron  Microscope  (SEM) 
micrographs.  A  substantial  AE  activity  was  characteristic  of  all  samples  tested  regardless  of 
the  character  of  their  macroscopic  brittleness  or  ductility  (Wong,  1990).  Thus,  by  inference,  a 
specimen  is  macroscopically  brittle  if  the  lateral  confinement  is  insuflident  to  keep  the  micro¬ 
crack  growth  stable.  Conversely,  if  the  microcrack  growth  is  stable  a  spedmen  is  macroscopically 
“ductile”. 

The  series  of  papers  hsted  in  the  above  paragraph  identified  two  major  mechanisms  of  microcrack 
nudeation  and  growth.  Hertzian  contact  stresses  between  caldte  cement  and  quartz  grairu 
(elastic  mismatch)  were  identified  as  the  major  micromechanical  process  leading  to  nurleation 
of  cracks  found  in  hydrostatically  compressed  candstone  samples.  The  contact  ctresses  were 
determined  by  modeling  the  porous  rock  as  a  randox^Jy  packed  assemblage  of  spherical  partides. 
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Since  the  contact  stress  concentrations  were  highly  localized,  the  initial  microcrack  length  was 
found  to  be  very  short  (ranging  from  0.004  to  0.3  mm).  The  final  crack  lengths,  radiating  from 
the  points  of  contact,  seemed  to  be  up  to  1  mm  long  (see  Fig.  6  in  Zhang,  et  al.,  1990a). 
The  second  micromechanical  mechanism  of  inelastic  rock  deformation  was  attributed  to  stress 
concentrations  at  the  perimeter  of  voids.  Sammis  and  Ashby  (1986)  have  developed  an  analytical 
model  for  the  latter  mechanism  based  on  fracture  mechanics  considerations,  which  sufficed  for 
formulation  of  relations  among  the  fracture  toughness,  the  porosity,  and  the  stress  field  (more 
specifically  the  ratio  between  the  confining  pressure  and  the  maximiim  compressive  stress). 

Additional  data  regarding  the  SEM  observations  of  microcrack  density  are  available  in  Wong,  et 
al.  (1992).  The  crack  densities  in  two  orthogonal  directions  were  determined  using  conventional 
stereological  methods  (Underwood,  1970,  Wong  and  Biegel,  1985).  In  hydrostatically-loaded 
Berea  sandstone  specimens,  the  number  of  microcracks  in  two  orthogonal  planes  was  almost 
identical  (with  differences  accounted  for  by  bedding);  however,  a  marked  increase  of  microcracks 
parallel  to  the  maximum  compressive  stress  emerged  when  a  differential  stress  was  added  to  the 
already  existing  hydrostatic  pressure.  Thus,  the  deformation  process,  reflecting  the  sequence  of 
different  modes  in  which  the  microstructure  changes,  strongly  depends  on  the  manner  in  which 
the  loads  are  applied.  In  other  words,  the  deformation  depends  on  whether  the  applied  loading 
is  proportional  or  not.  In  the  language  of  applied  mechanics,  the  deformation  processes  in  the 
rock  are  strongly  path  dependent. 

Precise  data  related  to  the  angular  distribution  of  cracks  are  even  more  difficult  to  find.  Hall- 
bauer,  et  al.  (1973)  carefully  examined  specimens  of  ar^aceous  quartzite  and  found  that  the 
majority  of  cracks  were  oriented  within  lO'’  to  the  compression  ajds*.  Microcracks  were  rather 
short  (0.1  to  1  mm  in  length),  intrsigranular,  and  confined  to  single  quartz  grains.  A  similar 
conclusion  was  arrived  at  by  Zheng,  et  aL  (1991)  who  investigated  compressive  stress-induced 
microcracks  in  limestone.  Most  microcracks  were  formed  by  bending  of  long  beam-like  grains. 
In  the  case  of  large  lateral  confining  stress  (34  MPa)  the  microcracks  were  reasonably  uniformly 
distributed  over  the  entire  volume  of  the  specimen  even  in  the  post-peak  regime.  At  a  lateral 
confining  stress  of  17  MPa,  localization  into  a  few  shear  bands  became  apparent.  At  zero  confin¬ 
ing  stress  the  damage  in  the  post-peak  regime  was  stron^y  localized  into  a  single  crack  near  the 
surface  of  the  specimen  parallel  to  the  compressive  axis.  The  rest  of  the  specimen  was  in  this  case 
almost  completely  undamaged.  According  to  the  data  in  Zheng,  et  aL  (1991)  the  average  angle 
subtended  by  microcrack  planes  and  the  compressive  axis  did  not  exceed  IS**.  The  microcrack 
density  (defined  as  the  number  of  cracks  per  unit  area)  was  found  to  increase,  and  the  microcrack 
length  decrease,  with  an  increase  of  confining  stress.  The  crack- width  to  crack-length  aspect 
ratio  was  typically  0.02.  Somewhat  older  data  for  sandstone  and  limestone  specimens  are  avsil- 
able  in  Swolfs  (1972),  Sangha,  et  al.  (1974),  Olsson  (1974),  and  Conrad  and  Friedman  (1976). 


*  See  page  vu  for  conversion  to  SI  units. 
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Application  of  acoustic  emissions  to  detect  the  location  and  intensity  of  microcracking  in  rocks 
has  been  suggested  in  the  recent  past  by  Holcomb  and  Costin  (1986),  Holcomb,  et  al.  (1990), 
Lockner  and  Byerlee  (1991),  Lockner,  et  al.  (1992),  and  others.  According  to  these  data  (see 
Lockner,  et  a/.,  1992)  A£  clustering  was  observed  in  Berea  sandstone  from  the  earliest  stages 
of  loading.  Initially  this  clustering  was  found  to  be  diffuse  suggesting  that  the  microcracks  can 
grow  only  through  Interaction  with  adjacent  defects.  The  level  of  localization  was  found  to  be 
directly  proportional  to  the  applied  compressive  force. 

In  summary,  the  available  data  set,  while  not  complete,  suffices  to  draw  a  number  of  impor¬ 
tant  conclusions  related  to  the  inelastic  deformation  of  sandstones.  At  low  temperatures  and 
stresses  below  the  crushing  level,  the  inelastic  deformation  of  sandstone  specimens  is  directly 
attributable  to  the  nucleation,  growth,  and  stability  of  microcracks.  The  specific  mechanism 
of  microcrack  nucleation  varies  with  microstructure.  The  growth  of  a  single  microcrack  (which 
does  not  interact  with  other  microdefects),  generated  from  a  pore,  appears  to  be  stable.  Thus, 
the  microcrack  interskction  problem  appears  to  be  instrumertal  in  micromechanical  modeling  of 
microcrack  growth.  The  distribution  of  microcrack  sizes,  shapes,  and  orientations  depends  on 
the  microstructure  (sizes  of  grains,  pores,  etc.)  and  the  stress  field.  The  brittle-ductile  transition 
and  the  final  failure  mode  depend  on  the  level  of  confinement. 

At  this  point  it  can  be  speculated  that  the  available  micromechanical  data  suffice  to  form  the 
basis  for  an  analytical  model  in  the  case  of  proportional  loading,  i.e.,  as  long  as  the  defor¬ 
mation  process  is  dependent  on  a  well  d.ocamented  sequence  of  dominating  micromechanisms. 
More  specifically,  during  proportional  loading  the  active  microcracks  will  likely  remun  active. 
Consequently,  changes  of  the  effective  macroparameters  (such  as  stiffiiess,  elastic  moduli,  and 
accumulated  damage)  will  change  gradually  and  in  orderly  fashion.  This  tjrpe  of  problem  is, 
for  simplidty,  best  handled  by  deformation-type  theories.  In  the  case  of  non-proportional  load¬ 
ings  characterized  by  rotation  of  the  principal  stress  directions,  the  available  microstructural 
experimental  data  must  be  appropriately  augmented  in  concert  with  analyses  and  macroscopic 
observations.  In  non-proportional  loading  a  group  of  microcracks  which  were  subjected  to  ten¬ 
sion  may  ultimately  find  themselves  in  hydrostatic  compression.  A  sudden  reversal  of  thrir  status 
firom  active  to  passive  will  result  in  a  discontinuous  change  of  the  effective  macroparameters. 
This  class  of  problems  can  be  handled  only  by  the  rate  theories. 

The  intricacies  of  the  brittle-ductile  transition  and  the  influence  of  the  rock  fabric  on  this  tran¬ 
sition  have  been  discussed  in  considerable  length  in  the  existing  literature  (see,  for  example, 
Paterson,  1978).  Analytical  models  of  this  transition  are  still  in  their  infancy.  TL's  is  not 
surprising  since  the  onset  of  the  brittle-ductile  transition  depends  on  the  microcrack  interac¬ 
tion.  The  extent  of  the  interaction  depends  on  the  distance  between  pairs  of  microdefects.  The 
largest  interaction  occurs  when  the  distance  between  two  defects  is  smallest.  In  other  words,  it 
is  the  smallest  distance  rather  than  the  average  distance  that  is  of  interest.  This  rather  simple 
conclusion  seems  to  have  been  totally  ignored  in  most  of  the  literature. 
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2.2.  CYLINDRICAL  PORE  MODEL 


In  porous  rocks  (such  as  sandstones  and  limestones),  microcrack  nucleation  is  commonly  at¬ 
tributed  to  the  local  tensile  stresses  along  the  pore  circumference.  The  circumferential  (hoop) 
normal  stresses  along  the  exterior  of  a  single,  circular  hole  of  radius  R,  embedded  in  a  plate  of 
homogeneous,  isotropic,  linear  elastic  material  (Fig.  2.1)  are 

(1  +  ^  -  (1  +  3^)co82tf]  ,  for  r  >  R  ,  (2.1) 

where  is  the  magnitude  of  the  macrostress  applied  in  the  direction  of  the  x-axis.  From  (2.1) 
the  extreme  values  of  the  stresses  at  the  pore  perimeter  are 


<Tg  =  3<r®  ,  for  ®  =  r  >  r  =  R 

A 

at  =  — <r®  ,  for  =  0  ,  r  ^  R  . 


(2.2) 


It  is  trivial  to  show  that  the  hoop  tensile  stresses  along  the  pore  drcumference  occur  within 
the  range  |  9  !<  30°.  However,  the  already  mentioned  observations  of  Hallbauer,  et  al.  (1973) 
on  argillaceous  quartzite  and  the  Zheng,  et  al.  (1991)  measurements  on  limestone  are  not 
in  very  good  agreement  with  the  results  predicted  by  this  simple  model.  According  to  thdr 
observations  there  are  no  apparent  cracks  with  orientations  larger  than  (beyond)  the  an^  of 
$  15°.  However,  a  crack  nucleated  in  a  plane  subtending  a  large  angle  to  the  compression 

axis  will  kink  and  grow  along  the  x-axis  (Fig.  2.1).  Moreover,  the  cracks  nucleated  at  9  s  0, 
t  —  R  will,  in  general,  become  unstable  before  the  stress  near  9  =  30°  becomes  large  enough  to 
propagate  cracks  from  small  notches.  Thus,  even  though  some  cracks  may  indeed  nucleate  from 
a  notch  at  an  angle  dose  to  9  =  30°,  the  average  crack  density  in  these  planes  will  be  minimal. 


Figure  2.1.  CyUndiical  pore  of  radius  R  under  remote  compressive  stress  <r°. 
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In  the  case  of  biaxial  compression 


(Tj;  =  -(T°  and  ay  =  -y;  <r®  ,  (2.3) 

where  ^  is  an  arbitrary  number,  and  the  tensile  stresses  vanish  for  <p>  1/3.  Thus,  according  to 
this  simple  criterion  the  microcracks  will  not  nucleate  at  the  perimeter  of  an  isolated  cylindrical 
pore  if  ^  >  1/3.  Introducing  the  average  (hydrostatic)  stress  p  and  the  deviatoric  (differential) 
stress  q 


P  = 


+  , 


and 


g  = 


the  local  tensile  stresses  wiU  vanish  if  the  ratio 


(2.4) 
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P 


l  +  9> 


<  1  . 


(2.5) 


The  equality  sign  in  (2.5)  provides  a  simpleminded  estimate  of  the  brittle-ductile  transition 
neglecting  microdefect  interaction,  other  micromechanical  mechanisms,  and  the  three- 
dimensional  character  of  the  phenomenon. 


It  is  also  of  interest  to  consider  the  stress  concentrations  in  the  proximity  of  two  interacting  cylin¬ 
drical  pores.  The  avidlable  data  in  Savin  (1961)  and  Paterson  (1978)  indicate  that  the  stress 
concentrations  are  decreased  by  the  presence  of  the  second  pore.  Detailed  rigorous  analyses  were 
reported  for  the  case  shown  in  Fig.  2.2  using  the  elasticity  model  formulated  by  Kouris  and 


Figure  2.2.  Two  Interacting  cylindrical  pores  under  remote  compressive  stress  o^. 
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Tschticida  (1991).  The  local  stresses  at  points  A,  B  and  C  are:  =  0.903ao  ,  =  -2.59cro  i 

and  Oy  =  O.SScTo  .  The  corresponding  values  for  a  single  pore  are  1.0  ,  —3.0  ,  and  1.0  ,  respec¬ 
tively.  However,  very  near  the  point  C,  i.e.  at  z  =  0.41,  the  matrix  is  subjected  to  biaxial 
tension  Og  =  O.OOoo  >Ad  Oy  =  0.36<ro  .  The  performed  calculations  indicate  that  the  second 
pore  has  a  shielding  effect,  and  will  seldom  if  ever  amplify  the  stresses  determined  for  a  single 
hole. 

Naturally,  the  perimeter  of  a  void  in  rock  material  is  never  smooth.  The  void  surface  is  typically 
irregular,  containing  small  notches  and  fissures.  Thus,  to  determine  the  stability  of  the  defor¬ 
mation  process  it  is  more  appropriate  to  consider  a  notched  void  as  shown  in  Fig.  2.3  (Kemeny 
and  Cook,  1991).  When  the  crack  length  a  is  much  smaller  than  the  pore  radius  R,  the  actual 
geometry  can  be  approximated  by  an  edge  crack.  This  edge  crack  is  subjected  to  the  local  hoop 
stress  computed  from  (2.2).  The  stress  intensity  factor  for  this  case  is  (Rooke  and  Cartwright, 
1976) 


Ki  =  1.12  (a,  -  Zoy)y/ra  ,  >  3ery  ,  a  «  R  .  (2.6) 

For  longer  cracks,  Sammis  and  Ashby  (1986)  suggested  an  approximate  expression  for  the  stress 
intensity  factor  in  the  form 


(11  .  (2-7) 

where  a  -  a/R.  The  normalized  strns  intensity  factor  <^y/itRlKc  »  where  Kc  is  the  critical 
stress  intensity  factor,  is  plotted  in  Fig.  2.4  vs.  the  aspect  ratio  a  s  a/R  ,  for  several  different 
confinement  ratios  ip  s=  .  In  all  cases  the  maximum  stress  intensity  factor  occurs  for  a  <R. 
After  the  stress  reaches  maximum,  the  crack  growth  becomes  unstable  (for  low  confinement)  in 
the  sense  that  the  crack  length  increases  for  decreasing  stress. 


Figure  2.S.  A  notched  void  of  radius  R  and  crack  length  a. 
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Sunmis  and  Aihby  (1986)  also  suggested  an  approximate  expression  for  the  case  of  interacting 
cracks.  According  to  this  suggestion  the  total  stress  intensity  factor  can  be  approximated  by 
adding  an  additional  term  to  (2.7) 

K',  =  —  {[1  -  ^(l  +  <«)’/lll  -  ^(l  +  k)V1  /(I  +  (2.8) 

T  TT  T 

reflecting  the  enhancement  attributable  to  the  direct  crack  interaction.  In  (2.8),  f  is  used  to 
denote  the  initial  porosity  of  the  rock. 


Figure  2.4.  Normalised  itrefs-istensity  factor  vs.  aspect  ratio  a/E,  for  several  different  con> 
finement  ratios 


2.3.  BRITTLE-DUCTILE  TRANSITION 

TJsbg  expressicms  (2.7)  and  (2.8)  the  Griffith’s  stability  condition  (Kanninen  and  Popdar,  1985) 
can  be  written  as 


Ki-\-K'jzzKc ,  (2.9) 

where  Kc  u  the  critical  stress  intensity  factor  considered  to  be  a  material  parameter.  The 
normalised  stress  needed  to  propagate  the  crack  can  now  be  written  from  the  three  above 
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expreuions  as  a  function  of  the  length  ratio  a  =  a/ R  and  the  degree  of  the  lateral  confinement 
^  >  0.  These  results  (taken  from  Wong,  1990)  are  reproduced  in  Fig.  2.5,  and  indicate  that 
confinement  as  low  as  ^  >  0.005  suffices  to  render  the  crack  stable.  This  conclusion,  however, 
contradicts  the  experimental  results  referenced  in  Wong  (1990)  indicating  that  the  Sammis  and 
Ashby  (1986)  model  overestimates  the  stabilizing  effect  of  the  lateral  confining  stresses,  i.e., 
that  it  proArides  a  low  estimate  of  the  magnitude  of  the  deviatoric  stresses  needed  to  prevent  the 
brittle  failure. 

Isida  and  Nemat-Nasser  (1987)  performed  rigorous  computations  which  indicated  that  the  ex¬ 
pression  (2.9)  for  the  stress-intensity  factor  overestimates  the  destabilizing  effect  of  the  crack 
interaction  (Wong,  1990).  Moreover,  the  two-dimensional  approximation  of  a  three-dimensional 
problem  further  exaggerates  the  influence  of  the  interaction.  Thus,  even  though  the  Sam¬ 
mis  and  Ashby  (1986)  model  exaggerates  the  effect  of  the  crack  interaction,  it  simultaneously 
overestimates  the  stabilizing  influence  of  the  lateral  confinement.  Thus,  this  model  leads  to  a 
contradictory  set  of  conclusions.  Even  though  their  interaction  is  exaggerated,  the  cracks  are 
still  predicted  to  grow  in  a  stable  mode  at  minute  levels  of  lateral  confinement.  This  apparent 
paradox  was  further  discussed  in  Wong  (1990).  The  influence  of  the  pore  shape  was  shown  to  be 
marginal  in  respect  to  the  determination  of  the  brittle-ductile  transition.  The  stress-intensity 
factors  for  a  crack  emanating  from  a  vertex  of  a  square  hde  (Hasebe  and  Ueda,  1980)  are  in¬ 
deed  larger,  but  not  enough  to  affect  the  stability  of  the  crack  growth.  The  expressions  for  a 
three-dimensional  case  of  spherical  cavities  (Sammis  and  Ashby,  1986)  do  not  change  the  basic 
conclusions  either.  This  prompted  Wong  (1990)  to  conclude  that  the  “Sammis  and  Ashby  (1986) 
model  is  not  appropriate  for  low-poronty  rocks”,  and  that  “compaction  mechanisms  (including 
pore  collapse  and  grain  rotation)  are  probably  operative”.  However,  it  semu  that  a  critical 
evaluation  of  the  pore  model  may  frimish  a  new  insight  into  the  phenomenon  and  provide  a 
better  estimate  of  the  brittle-ductile  transition. 

The  Isida  and  Nemat-Nassa  (1987)  computations  were  performed  for  a  perfect,  doubly  periodic 
pattern  of  two  cracks  emanating  from  a  circular  hole.  Using  these  results  to  explain  experimental 
observatioiu  regarding  the  onset  of  unstable  crack  growth  implies  self-similar  growth  of  defects 
throughout  the  deformation  process.  However,  as  shown  in  Bazant,  et  aL  (1989)  and  Bazant 
and  Cedolin  (1991)  self-similar  growth  of  defects  represents  a  thermodynamically  unstable  path. 
This  is  obvious  on  purdy  physical  grounds  as  weU.  Assuming  all  distances  between  adjacent 
defects  to  be  equal  (as  implied  by  the  so-called  cell  method  based  on  doubly  periodic  patterns) 
the  minimum  distance  between  two  defects  is  obviously  maximized.  However,  the  interaction  of 
cracks  and  the  formation  of  a  cluster  of  cracks  depends  almost  entirely  on  the  TninitnnTn  (rather 
than  average)  distance. 
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Figure  2.5.  Nonnalized  stress-inteiuity  factor  ti.  aspect  ratio  a/R,  for  several  different  values 
of  porosity  /  and  two  confinement  ratios  tp. 
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The  already  mentioned  acoustic  emission  data  on  sandstones  (Lockner,  et  al.,  1992)  provides 
strong  evidence  that  the  A£  signal  clustering  characterizes  the  deformation  process  from  its  ear¬ 
liest  stages.  Most  of  the  signals  are  strongly  clustered  indicating  strong  microcrack  interaction. 
This  experimental  evidence  proves  the  existence  of  a  basic  flaw  in  the  predictions  of  the  cell 
models  which  imply  preservation  of  the  periodicity  of  the  defect  pattern.  The  assumption  that 
the  defect  periodicity  (microstructural  order)  perseveres  into  the  latter  part  of  the  deformation 
process  dominated  by  the  microdefect  interaction,  substantially  overestimates  the  ductility  and 
strength  and  underestimates  the  onset  of  brittle  instability.  This,  in  fact,  renders  the  cell  mod¬ 
el  inadequate  in  qualitative  and  qus  atitative  senses  alike.  It  seems  reasonable  to  acknowledge 
the  dominsmt  role  of  the  clustering,  i.e.,  the  spatial  disorder  of  the  microstructure.  Since  the 
clustering  at  low  density  of  defects  depends  on  the  smallest  distance  between  adjacent  cracks, 
it  is  important  to  investigate  the  effect  of  the  statistical  distribution  of  defect  distances  on  the 
deformation  pattern. 

It  is  important  to  emphasize  that  the  available  experimental  data  clearly  indicate  that  the  onset 
of  the  brittle-ductile  transition  is  a  statistical  event.  The  emergence  of  the  zone  of  large  micro¬ 
crack  density  depends  not  only  on  the  porosity  (first  statistical  moment  of  the  void  distribution), 
but  in  a  much  more  substantial  manner  on  the  higher,  perhaps  extreme,  statistical  momenta 
defining  extremes  of  the  pore  spacing.  Thus,  a  study  «iTwnaT  to  one  suggested  by  Horii  and 
Nemat-Nasser  (1986)  for  compact  rocks  may  not  provide  a  reliable  estimate.  In  fact,  it  seems 
that  every  deterministic  modd  based  on  the  average  spacing  of  pores  will  underestimate  the 
onset  of  the  unstable  crack  regime. 
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S.O  PHENOMENOLOGICAL  DAMAGE  MODELS  FOR 
BRITTLE  DEFORMATION  OF  ROCKS 


3.1.  CRACK  DENSITY  DISTRIBUTION 

Consider  n  porous  rock  specimen  containing  a  certain  distribution  of  microcracks  accumulated 
during  a  specific  loading  program  from  some  initial  state.  Two  questions  of  paramount  impor¬ 
tance  are:  a)  wbat  are  the  damage  variables  that  adequately  represent  the  degraded  state,  and 
b)  how  does  the  elastic  stiffness  of  a  damaged  rock  specimen  depend  on  the  introduced  damage 
variables?  Regarding  the  first  question,  various  damage  variables  have  been  introduced  in  the 
literature.  Most  of  these  models  are  of  limited  validity.  If  the  current  crack  pattern  in  the  rep¬ 
resentative  volume  element  is  such  that  cracks  are  uniformly  distributed  in  all  planes,  regardless 
of  their  orientation,  a  scalar  damage  variable  should  be  a  natural  choice.  The  corresponding 
distribution  of  damage  is  referred  to  as  isotropic.  If  cracks  are  nonuniformly  distributed  over 
differently  oriented  planes,  damage  distribution  and  correspondingly  the  material  response  are 
anisotropic.  A  distribution  function  p(n)  (defined  on  a  unit  sphere)  can  be  introduced  to  define 
the  directional  dependence  of  the  crack  density.  This  function  can  be  expanded  in  a  Fourier- 
type  series  of  certain  families  of  spherical  functions  (Kanatani,  1984,  Onat  and  Leckie,  1988), 
containing  dyadic  products  of  the  unit  vector  and  the  Kronecker  delta  tensor.  In  addition  to  the 
scalar  (isotropic)  term,  the  second-,  fourth-,  and  higher  even-order  sym*netric  tensors  appear  in 
this  representation.  Therefore,  the  accurate  representation  of  a  complicated,  highly  anisotropic 
orientation  of  damage  by  introdudng  some  average  tensor  measure  of  damage  is  a  difficult  task 
which  in  general  requires  introduction  of  second-,  fourth-,  and  possibly  even  higher-order  ten¬ 
sors  to  represent  the  state  of  damage.  An  altunative  description  of  damage  anisotropy  involving 
stereological  measurements  and  using  a  geometric  probability  approach  is  discussed  by  Wong 
(1985). 

In  a  general  case  of  loading  of  initially  anisotropic  rocks,  the  planes  containing  extreme  densities 
of  damage  are  not  mutually  perpendicular.  Consequently,  both  the  damage  itself  and  its  effect 
on  the  material’s  effective  sti&ess  are  anisotropic.  However,  in  the  case  of  initially  isotropic 
brittle  rocks  subjected  to  proportional  loading,  the  density  of  damage  is  maximum  in  the  plane 
perpendicular  to  the  largest  prindpal  tensile  stress  and  minimum  in  the  plane  normal  to  the 
wnifiimiiTn  principal  stress.  To  study  this  case  of  damage  distribution,  it  seems  reasonable  to 
approximate  its  density  distribution  by  an  oval  (Fig.  3.1a).  This  type  of  damage  distribution 
can  then  be  represented  by  a  second-order  tensor.  H  pij  denotes  the  components  of  the  second- 
order  crack  density  tensor  p,  the  density  of  cracks  embedded  in  the  planes  with  a  normal  n  is 
given  by 


p(n)  =  pki  . 


(3.1) 
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(C) 


Figure  S.l.  (a)  Crack  diitribution  within  orthogonal  crack  famil*^  modeled  by  a  continnons 
crack  distribution  of  oval  shape;  (b)  If  the  crack  densities  are  the  same,  the  oval  shape  becomes 
spherical,  giving  rise  to  an  isotropic  damage  distribution;  (c)  Crack  distribution  corresponding 
to  inclined  crack  families. 
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depicted  by  the  oval  shape  of  Fig.  3.1a.  Multiplying  (3.1)  with  n,nj  and  integrating  over  all 
directions  spanning  the  entire  solid  angle  fl  =  4ir  ,  and  using 


/  ninjUktii  dfl  =  — 


it  follows  that 


8ir  1  / 

•  (3.2) 

Above,  S^j  denotes  the  Kronecker  delta.  The  contraction  t  =  y  in  (3.2)  defines  the  first  invariant 
of  the  crack  density  tensor 


Pkk  = 


Substituting  (3.3)  into  (3.2),  the  crack  density  tensor  can  be  expressed  as 


Pij  —  ~  5F0  ^ij)  • 


It  (3.4) 


if  the  density  of  aU  cracks  within  a  unit  volume,  and 


(3.3) 


(3.4) 


(3-6) 


^  (3.6) 

is  a  second-order  tensor  which  shall  be  referred  to  as  a  damage  tensor.  The  corresponding  damage 
distribution  is  orthotropic.  If  the  distribution  is  isotropic,  so  that  p(u)  s  Po/4t  =  constant , 
(3.6)  reduces  to  Vij  =  \pofij  >  while  the  crack  density  tensor  (3.4)  becomes  pij  =  ^Po^ij  • 
The  oval  shape  of  Fig.  3.1a  then  transforms  to  a  spherical  shape  (Fig.  3.1b). 

In  two-dimensional  problems  one  has 


2v 


i 


SO  that,  in  place  of  (3.4),  the  crack  density  tensor  is  ^ven  by 
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(3.4') 


Pii  =  “(®u  “  aPp  • 

Higher-order  tensors  are  needed  to  accurately  represent  the  general  cases  of  crack  distributions, 
especially  when  the  damage  is  not  orthotropic.  For  example,  the  fourth-order  crack  density 
tensor  Pijui  describes  the  crack  density  corresponding  to  planes  with  a  normal  n  as 


p(n)  =  pijki  nin^ntn,  .  (3.7) 

The  second-order  damage  tensor  is  still  defined  by  (3.6),  while  the  fourth-order  damage  tensor 
is  defined  by 


=  /  p(n)ni»i»fc»/  dn  .  (3.8) 

•'4* 

The  relationship  between  the  fourth-order  crack  density  tensor  p  and  the  second-  and  fourth- 
order  damage  tensors  2>  and  T)  can  be  derived  by  a  similar  procedure  as  (3.4).  It  follows  that 

Pm  =  +  ^pp »««) .  (3-8') 

where: 


Aijki  —  +  SikDji  -1-  SiiVjk  +  +  iji^ik) 

Bijkl  =  +  fik^jl  +  dijfjfc)  • 

In  two  dimensions,  expression  (3.8')  is  replaced  by 

Pijki  =  ^(Pijkl  -  Bijkl)  •  (3.8") 

3.1.1.  Second-  And  Fourth-Order  Approximations 

To  illustrate  application  of  the  results  from  Section  3.1,  consider  first  the  case  of  two  orthogonal 
crack  families,  each  having  the  crack  density  ^/2.  In  a  two-dimensional  analysis  this  crack 
distribution  can  be  represented  by 

^*)  =  +  <(»  -  |)  +  «{*  -  x)  +  <(*  -  ^)1  .  (a) 
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where  B  defines  aui  arbitrary  direction  through  the  point  under  consideration,  and  6  denotes  the 
Dirac  delta  function.  Substituting  (a)  into  the  expression  for  the  second-order  damage  tensor 


Vij  =  /  p(n)n.nj  dil  ,  (3.6') 

•'2ir 

it  follows  that  Vij  =  \po6ij  .  Hence,  from  (3.4')  the  crack  density  tensor  is  pij  ~  ^po^ij  , 
so  that  within  the  second-order  approximation  (3.1),  the  distribution  (a)  is  replaced  by  the 
continuous  distribution 


p{e)  =  PijTunj  =  —po  ,  (6) 

i.e.,  the  circle  of  radius  Pq/2t  .  Therefore,  two  orthogonal  crack  families  with  the  same  crack 
densities  are  replaced  by  an  isotropic  homogeneous  crack  distribution. 


A  more  accurate  approximation  is  obtained  using  the  fourth-order  damage  tensor.  Substituting 
(a)  into  (3.8),  it  follows  that  Dim  =  D2223  =  po/^  ,  while  all  others  components  are  equal  to 
zero.  The  tensor  A  has  the  components  Xmi  =  .A2222  =  Pol^  ,  a^d  An22  =  >^2211  =  •A1212  = 
•A2121  =  ^1221  =  *^2112  =  Po/6  I  as  the  only  non-zero  components.  Substitution  into  (3.8")  and 


Figure  S.2.  Second-  and  fourth-order  continuous  approsdmations  of  two  orthogonal  crack 
families  with  the  same  crack  densities. 
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(3.7),  therefore,  gives  the  continuous  fourth-order  approximation 


P{^)  =  Pijkininjnknt  =  (1  +  2cos4d)  .  (c) 

The  second-  and  fourth-order  continuous  approximations  (b)  and  (c)  are  plotted  in  Fig.  3.2. 
A  remarkable  feature  of  the  fourth-order  approximation  in  this  case  is  that,  in  addition  to 
dominating  regions  of  positive  crack  density,  two  orthogonal  regions  of  negative  crack  density 
also  occur,  at  45°  relative  to  the  positive  regions.  This  can  even  happen  within  the  second-order 
approximation,  as  will  be  discussed  in  the  next  subsection. 

3.1.2.  Negative  Crack  Density— Antieraeks 

The  observed  feature  of  negative  crack  density  often  occurs  when  the  actual  crack  distribution 
is  approximated  by  a  continuous  distribution  corresponding  to  second-,  fourth-  or  higher-order 
damage  tensors.  To  illustrate  this,  consider  a  single  family  of  cracks  whose  normal  is  n  =  {0, 1} , 
written  in  a  symmetrical  form  as 


p(*)=fW*-f)+<(*-^)J-  (<0 

The  continuous  approximation  of  this  distribution,  corresponding  to  the  use  of  the  second-order 
damage  tensor,  is 


=  ^Po(l  -  2  cof2^)  .  (e) 

Derivation  of  the  above  expression  is  identical  to  one  already  explained  in  Section  3.1.1.  The 
graphical  depiction  of  the  second-order  tensor  approxiination  is  shown  in  Fig.  3.3.  The  interest¬ 
ing  feature  of  this  continuous  distribution  is  emergence  of  the  negative  crack  density  over  a  part 
of  the  range.  This  means  that  in  the  corresponding  regions,  actual  cracks  are  replaced  by  rigid 
lameDa  (the  opposite  of  a  crack)  which  are  referred  to  as  negative  cracks  or  antieraeks  (Dundurs 
and  Markenscoff,  1989).  The  emergence  of  negative  crack  densities  as  a  result  of  approximating 
discontinuous  distributions  of  cracks  by  continuous  distributions  represented  by  tensors,  should 
have  been  expected.  A  tensorial  appnudmation  (e)  of  a  delta  function  implies  the  existence  of 
damage  at  angles  other  than  0  =  ir/2  and  $  s  3ir/2.  Consequently,  negative  crack  densities 
must  be  present  to  balance  this  nonexisting  damage. 

The  regions  with  negative  crack  density  and  corresponding  anticrack  distributions,  also  occur 
when  the  fourth-order  approximation  is  used,  bi  this  case  one  has 
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(/) 


p{0)  =  ;^Po(l  -  2  cos20  +  2  cos4^)  . 

The  plot  of  this  density  distribution,  with  its  regions  of  negative  crack  density,  is  shown  in  Fig. 
3.3  (dashed  curve).  The  actual  crack  distribution  is  approximated  much  better  than  in  the  case 
of  the  second-order  tensor. 


Figure  S.S.  Second-  and  fourth-order  continuous  approximations  of  a  single  crack  family. 
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Of  course,  iu  many  cases  the  approximate  continuous  crack  distribution  does  not  contain  regions 
of  negative  crack  density.  For  example,  for  two  orthogonal  crack  families  with  crack  densities  pi 
and  P2  —  Po  —  Pi  I  where  po  is  the  total  crack  density,  the  occurrence  of  regions  with  negative 
crack  densities  depends  on  the  ratio  pi/p^  .  If  pi  =  po/8  and  pi  =  Tpo/d  ,  the  density 
distribution  according  to  the  second-order  approximation  is 


P{^)  =  “  ^cos2d)  . 

Thus,  the  negative  crack  density  regions  occur  ag^,  as  shown  in  Fig.  3.4a.  However,  if 
Pi  =  pol2  and  p}  =  2po/3  ,  the  crack  density  distribution  is 


showing  no  regions  of  negative  crack  density  (Fig.  3.4b).  The  tramsition  case  is  pi  =  po/4  and 
P2  =  3po/4  ,  for  which 


P{9)  =  ^^(1  -  , 


(see  Fig.  3.4c). 


(a)  (b) 


(C) 


Figure  S.4.  Continuous  crack  ^tributions  corresponding  to  two  orthc^onal  crack  families  of 
different  crack  density  ratios. 
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3.1.3.  Rose  Diagram 


Tlie  same  procedure  applies  in  approximating  an  arbitrary,  experimentally  measured  crack  dis¬ 
tribution  by  a  continuous  distribution,  based  on  the  even-order  tensor  measures.  Consider  for 
example  a  “rose  diagram”  deduced  from  the  actual  measurements  in  tests  performed  by  Hall- 
bauer,  et  aJ.  (1973),  Fig.  3.5a.  Using  the  second-order  tensor  approximation,  the  components 
of  the  corresponding  damage  tensor  are  calculated  from  (3.6')  with  n  =  {cosd.sin^}  to  be  : 
Pii  =  0.061po  >  ^22  =  0.939/£>o  aod  ^12  =  O.lllpo  •  Consequently,  substituting  these  results 
into  (3.4')  the  components  of  the  crack  density  tensor  are  computed  to  be:  pn  =  -^0.75po  , 
P22  =  ^2.75po  I  and  pu  =  2^0.44po  •  Therefore,  the  approximate  continuous  crack  distribution 
is  given  by 


i.e. 


I 


P(^)  —  ;^Po(— 0-75  cos^ff  -f  2.75  sin*^  -f-  0.44  sin0cos0)  , 
2ir 


p(ff)  =  — /)o(l  +  0*22  sin2d  -  1.75  cos2^)  . 


Figure  3.S.  (a)  A  “rose  diagram”  deduced  from  actual  test  measurements.  The  representative 
crack  densities  are:  pi  =  =  1.8po/2x  ,  p2  =  14Apol2v  ,  ps  =  10.8po/2T  ,  and  p4  = 

7.2piil2it  ,  where  po  is  the  total  crack  density;  (b)  Continuous  second-order  approximation  for 
crack  density  distribution  corresponding  to  the  above  “rose  diagram”. 
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This  approxiinate  crack  distribution  is  plotted  in  Fig.  3.5b.  In  addition  to  the  dominating 
regions  of  positive  crack  density,  the  regions  of  negative  crack  density  are  still  present.  The 
scalar  (isotropic)  approximation  is  a  circle  of  radius  po/2ir  .  This  phenomenon,  as  yet  not 
discussed  in  the  literature,  deserves  careful  consideration.  This  is  especially  true  since  negative 
crack  densities  actually  do  occur  in  analyzing  situations  such  as  those  studied  by  Hallbauer,  et 
al.  (1973). 

3.2.  REPRESENTATION  OF  DAMAGE 

To  Illustrate  a  class  of  relatively  simple  representations  of  damage,  consider  an  initially  isotropic, 
homogeneous  and  elastic  matrix  subjected  to  proportional  loading.  According  to  the  experimen¬ 
tal  data  on  rocks  (Hallbauer,  et  al.,  1973,  or  Zheng,  et  al.,  1991 )  the  planes  of  extreme  microcrack 
density  coincide  with  the  principal  planes  of  the  stress  tensor.  Consequently,  the  specimen  is 
orthotropic  on  the  macroscopic  scale  suggesting  that  the  microcrack  distribution  can  be  repre¬ 
sented  by  the  second-order  damage  tensor  (Vakulenko  and  Kachanov,  1971,  Kachanov,  1980, 
etc.) 


D  =  pnN  +  pnM  -I-  pkK  ,  (3.9) 

where  N  =  n®ii,  Msm®m  and  K  «=  k  ®  k  are  the  dyadic  products,  and  n,  m  and 
k  axe  the  unit  normals  to  three  crack  families.  The  scalar  quantities  pn,  Pm  and  pk  are  the 
corresponding  crack  densities. 

E^erimental  determination  of  crack  densities  involves  a  micro-to-macro  transition.  This  implies 
the  existence  of  a  rq>resentative  volume  element  centered  at  a  point,  which  contains  a  statistically 
valid  sample  of  microcracks  influendng  the  state  at  the  considered  pdnt.  Thus,  if  mthin  a 
material  element  of  unit  volume  there  are  N  parallel  flat  cracks  with  the  unit  normal  n  and  an 
average  characteristic  length  a,  the  nondimenidonal  crack  density  is  defined  as  pn  =  Ne?  .  The 
crack  length  a  is  calculated  as  an  average  characteristic  length  of  all  cracks  within  a  representative 
volume  element,  imbedded  on  planes  with  the  normal  n. 

The  minimmn  prescription  for  representation  (3.9),  based  in  equal  parts  on  intuition  and  ex¬ 
perimental  data,  is  to  measure  microcrack  densities  p^,  p^  and  pk  in  three  principal  planes. 
Naturally,  there  is  no  assurance  that  the  estimate  of  microcrack  density  in  an  arbitrary  plane, 
interpolated  using  the  representation  (3.9),  wiU  coincide  with  the  actual  density.  Indeed,  the 
representation  of  the  damage  distribution  by  the  second-order  tensor  (3.9)  amounts  to  assum¬ 
ing  a  continuous  crack  distribution  of  the  oval  type  shown  in  Fig.  3.1a,  with  appropriately 
adjusted  crack  densities  (preserving  as  constant  the  total  number  of  cracks  in  a  unit  volume, 
i.e.  Pn  +  Pm  +  ^  .  gi^ea  ky  (3.5)).  If  p„  =  p,„  =  p*  =  po  ,  the  oval  shape  becomes  a 

sphere  (of  radius  3po/4x).  Therefore,  within  a  description  by  the  second-order  damage  tensor. 
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the  three  orthogonal  crack  familiec  with  the  same  crack  densities  give  rise  to  an  isotropic  damage 
distribution  (Fig.  3.1b). 

A  more  ambitious  alternative  consists  of  measuring  microcrack  densities  p'  in  a  large  number  of 
planes  with  different  orientations  n',  in  order  to  form  a  “rose”  histogram  (Kanatani,  1984).  In 
this  case  the  second-order  representation  (3.9),  in  which  pn,  Pm  and  pk  are  the  principal  values, 
and  n,  m  and  k  are  the  principal  directions  of  the  second-order  tensor  p'n'0n'  ,  may  or  may 

not  be  a  satisfactory  approximation  of  the  measured  data.  Geometrically,  the  approximation 
(3.9)  again  amounts  to  replacing  a  more  complicated  distribution  (such  as  the  one  shown  in  Fig. 
3.1c)  by  an  oval  (“orthotropic”)  distribution. 

Since  N  -I-  M  -f  K  =  1  ,  (3.9)  can  be  rewritten  as 

®  =  {pn  —  ^i:)N  4-  {pm  —  P/t)M  -f  />fcl  ,  (3.10) 

where  1  denotes  the  second-order  unit  tensor.  Two  special  cases  of  (3.10)  are  of  interest.  If 
Pk  —  Pm^  Pn  (transversely  isotropic  case),  the  damage  tmsor  (3.10)  simplihes  to 


l>  =  (p«-/»m)N  +  p«l  .  (3.11) 

In  the  Isotropic  case,  p^  —  pm  —  Pk  —  P  t  vid  (3.11)  reduces  to  a  spherical  tensor  D  s  pi  , 
describing  an  isotropic  damage  distribution. 

3.3.  ELASTIC  STIFFNESS  TENSOR 

In  the  considered  brittle  deformation  processes,  the  state  of  the  material  is  locally  defined  by 
two  state  variables:  elastic  strain  £«  and  the  damage  tensor  V.  The  second  of  these  variables  is 
often  referred  to  as  the  internal  or  hidden  variable.  The  elastic  free  energy  iS|  therefore,  = 
fft«(E«,l>)  .  This  e3q>ression  must  be  invariant  with  respect  to  any  change  of  reference  frame. 
Since  under  the  change  of  frame  defined  by  an  orthogonal  transformation  Q,  the  damage  tensor 
V  becomes  Q2>Q^  (superscript  T  denotes  the  transpose),  the  invariance  condition  requh-es 
the  elastic  free  energy  to  be  an  isotropic  function  of  both  dastic  strain  Ee  and  damage  V. 
According  to  the  well  known  invariance  theorem  (Spencer,  1971),  an  isotropic  scalar  function  of 
two  symmetric  second-order  tensor  variables  can  be  represented  as  a  polynomial  of  its  irreducible 
integrity  basis: 


(E.  :  1),  (E.  :  Ee),  (E.’  :  Ee),  {V :  1),  (P :  P),  (P"  :  P) 
(Ee  :  P),  (Ee  :  P"),  (E**  :  P),  (Ee*  :  P*) . 
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In  the  cue  of  brittle  deformation  processes,  the  elutic  strains  are  typically  infinitesimal.  More¬ 
over,  it  will  be  usumed  that  the  stress  is  equal  to  sero  when  the  elutic  strain  vanishes.  Thus, 
the  elutic  free  energy  is  a  quadratic  function  of  the  strain  components 

•  Ee)  +  %(£*  :  1)^  +  »b(Ee  t  l)(Ee  ^ 

+  qs(E.*  :  D^)  +  tte(E,  :  Vf  +  t^{Ee  : 

In  (3.12),  tji  (t  =  1,2, 3,4)  are  the  constants  or,  more  generally,  the  scalar  functions  of  the 
invariants  of  X>.  Using  (3.12)  and  the  expression  tr  =  djptl^^e  ■  the  stress-strain  relationship 
acquires  the  following  form 


tr  —  2f;iEe  +  2i^(Ee  :  1)1  +  fti[(Ee  :  D)!  +  (E,  :  1)P]  +  *74(Ee  ■  P  +  D  •  E*) 

+  i75(Ee  •  P*  +  •  Ee)  +  2%(E,  :  P)P  +  2q7(Ee  :  P*)P^  . 

Above,  (•)  denotes  the  inner  tensor  product.  Expression  (3.13)  is  cleuly  lineu  in  strain,  admit¬ 
ting  the  familiar  form 


tr  =  :  Ee  .  (3.14) 

The  Cartesian  components  of  the  instantaneous  elutic  stifiEness  tensor  £«,  written  in  the  sym- 
metridxed  form,  axe 


^ijkl  —  +  SiiSjic)  +  2rf2  SijSkl  +  Vs{^ij‘f^kl  +  ^ij^kl) 

f  (S.16) 

+  2  +  ^il^jq^qk  +  ^if^fl^jk) 

+  2rie  VijVici  +  2iyr  t>ipVpjVifqVgt  • 

In  (3.15),  Sij  denotu  the  Kronecker  delta,  i.e.,  the  components  of  the  second-order  unit  tensor 
1.  Repeated  indices  indicate  summation.  Expression  (3.15)  possesses  the  obvious  symmetries 
^ijkt  =  ^jikl  —  ^ijlk  »  u  well  u  the  self-adjoint  symmetry  ^ijkl  s  .  Hence,  the  elastic 
strain  energy  (3.12)  can  be  written  u 

^.  =  i£e(l>):(Ee®Ee).  (3.16) 

Conuder  the  case  ^hmi  the  damage  tensor  is  defined  by  (3.11).  Since  NikNkj  =  Nij  =  rtiUj 
and  Nkk  =  nk’^k  =  1  ,  it  follows  that: 


SijT^kl  +  T>ijSkl  —  2pm  lijkl  +  i^ijkl  +  ^jkt) 
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iikT^jl  +  T^ik^jl  +  iuT^jk  +  ^uSjk  =  ^Pm  l}jkl  +  4Ap  /,*•*/ 


iikT>i<p9l  +  T^iq^qk^jl  +  SuVj^V^k  +  T^iq^ql^jk  =  ^Pm  ijjkl  +  ^^Pm  +  Ap)Ap  J,®  */ 
T>ijVkl  =  pin.  lijkl  +  Pm^P  {lijki  +  ^ijkl)  +  (^P)^  Ifjkl 

‘^ifJ>pj^kq^ql  =  Pm  ^ijkl  +  Pm(2p»n  +  Ap)Ap  {liju  +  lijkl)  +  (2pm  +  Ap)*(Ap)’  ifjh,  , 

where  Ap  =  Pn  —  Pm  •  The  tensors: 

-r?j«  =  +  SiiSjk)  ,  =  SijStt 

^ijkl  ~  ^ij^k^l  I  ^jkl  ~  (3*17) 

^jki  =  +  Siinjrik  +  mniSjk)  ,  -^ik'  ~  »!«;«*»/ 

combining  dyadic  (tensor)  products  of  the  Kronecker  delta  and  a  unit  vector,  form  the  integrity 
basis  for  the  fourth-order  tensors  that  are  symmetric  with  respect  to  the  first  and  second  p^ 
of  indices  (Kunm,  1983).  The  finear  tensor  space  spanned  by  this  basis  is  closed  with  respect 
to  the  trace  product,  forming  an  algebra.  For  example,  it  can  be  shown  that  P  :  P  =  P  , 
:  P  =  31^  ,  I’  :  P  =  and  I*  :  P  =  .  Hence,  the  stiffness  tensor  (3.15)  can  be  written 


A  =  2(in  +  V4pm  +  nspm)  I'  +  2(i7J  +  1j3Pm  +  VtPm  +  ^Pm)  I* 

+  [*J3  +  2»^Pm  +  27fTpln{2pm  +  Ap)] Ap  (I®  +  I*  )  (3.18) 

+  2(74  +  7s(2pm  +  Ap)]Ap  I®  +  2Ii}6  +  «77(2pm  +  Ap)*](Ap)’  I'  . 


3.3.1.  Jsatropie  Damage 

Consider  first  the  special  case  of  isotropic  damage.  In  this  case  the  damage  is  folly  defined 
by  a  single  scalar  variable,  i.e.,  Vij  =  V  %  ,  where  V  is  the  crack  density  in  each  plane. 
Geometrically,  the  oval  distribution  of  Fig.  3.1a  reduces  to  a  spherical  distribution.  Since  the 
crack  density  does  not  depend  on  the  orientation  of  the  plane,  p^  —  Pm  —  P  Ap  =  0  ,  and 
the  sti&ess  tensor  (3.18)  reduces  to  a  simple  form 
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(3.19) 


A  =  2(t7i  +  +  VsP^)  I*  +  2(^^2  +  rj3P  +  VeP^  +  rfjp*)  V  . 

Neglectmg  the  term  proportional  to  p*,  i.e.  taking  f77  =  0  ,  and  introducing  the  notations: 

rh  =Vs  =  Po  ,  Vi  =  -2mo 

wo  _  .  (3-20) 

tj2  =  »J6  =  ■*0/*  »  V3  — 

where  Ao  and  po  are  the  Lame  elasticity  constants  of  the  undamaged  material,  (3.19)  becomes 

£e  =  (l-  P?  {2pc  I'  +  Ao  I*)  =  (1  -  py  .  (3.21) 

In  (3.21),  £t^  =  2po  +  Aq  P  is  the  stiffness  of  the  virgin  material.  Defining  <<;  =  (1  -  p)^  as 
the  new  damage  variable,  (3.21)  takes  the  form  of  the  familiar  representation  of  the  isotropically 
degraded  elastic  stiffiiess,  £e=u  ,  originally  proposed  by  Kachanov  (1958)  and  considered 
at  length  by  Lemaitre  (1987,  1992). 

In  this  case  the  damage  distribution  is  fully  defined  by  a  single  parameter  p.  This  parameter  can 
be  readily  measured  by  comparing  slopes  of  the  unloading  segments  of  the  stress-strain  curves. 

3.3.2.  Transversely  Isotropic  Damage 

If  the  crack  distribution  is  not  isotropic,  but  approximated  by  a  second-order  damage  tensor 
1>,  the  parameters  17, ■  defined  by  (3.20)  should  be  accordingly  adjusted,  in  concert  with  the 
experimental  data  for  the  particular  material  sad  crack  arrangement.  In  the  sequel,  however, 
the  material  response  wiU  be  approximated  by  the  stiflhess  tensor  (3.18)  and  parameters  Vi  pven 
by  (3.20).  Hence 


r.  =  (1  -  Pm?  {2li0  I*  +  Ao  V)  -  (1  -  Pm)Ap  Ao  (I’  +  I^) 

(3.22) 

-  [2(1  -  p,„)  -  ApjAp  2aio  I' +  (Ap)»  Ao  I* 

defines  the  elastic  stifihess  tensor  of  the  considered  transversely  isotropic  material,  whose  axis 
of  rotational  symmetry  is  parallel  to  the  direction  n.  In  this  case  two  parameters  define  the 
damage  distribution:  pm  and  Ap  s  pn  "  />m  • 

H  the  crack  densities  are  sufildently  small,  rendering  the  terms  contuning  the  square  of  the 
crack  densities  negligible,  (3.22)  reduces  to 

£e  =  i:.®  -  2pm  (2/10  +  Ao  I’)  -  Ap  Ao  (I®  + 1^)  -  4Ap  /lo  I*  .  (3.23) 

Finally,  if  all  cracks  are  embedded  in  planes  parallel  to  n  ( pn  =  0  ,  Ap  =  — pm  )?  (3.23)  becomes 
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A  =  -/>m  [  4/io  I‘  +2Ao  I*  -  Ao  (V  +  V)-4  nol^]. 


(3.24) 


3.4.  ELASTIC  COMPLIANCE  TENSOR 

To  invert  the  elastic  stiffness  tensor,  expression  (3.22)  must  first  be  rewritten  as 

6 

/:,  =  5;  o.  r  ,  (3.25) 

«=i 

where  the  scalar  parameters: 

=  2/io(l  -  pmY  ,  =  Ao(l  -  pm)^ 

03  =  — Ao  (1  —  Pm)I^P  I  =  -Ao  (1  -  pm)^P  (3.26) 

as  =  -Ipa  [2(1  -  pm)  -  ApjAp  o®  =  Ao  (Ap)^ 

depend  only  on  the  material  properties  of  the  undamaged  material  (Ao,Po)  and  the  densities 
of  the  already  accumulated  damage.  To  invert  the  expression  (3.25),  i.e.,  to  derive  the  elastic 
compliance  tensor  M.,  s  ,  it  is  convenient  to  change  the  basis  I*  into  J*  (t  s  1,2,...  ,6) 


through  the  linear  transformation: 

J^si(l’-I»-I<  +  3I«), 

J®  =  i(-2I»-I"  +  3I«), 

J*  =  2(I*-I®), 

J*  =  i(_l2+I3  +  I4+I«) 

J<  =  |(-2I’  +  I<+I«) 

J«  =  1(2I*-I*+I®  +  I^-4I*+I®)  . 

(3.27) 

The  elastic  sti&ess  (3.25)  then  becomes 


6 

r.  =  j;  bi  J<  .  (3.28) 

«sl 

The  new  parameters  fr,-  in  (3.28)  are  related  to  the  parameters  Oj  in  (3.26)  by: 

=  5(2ai  +  3a2+2o3  +  Os+a6),  hj  =  l(-02  +  2a3  +  aj  +  a*) 

63  =  -|(a2  +  03)  ,  ^  ^  (3.29) 

1 

hi  S=  2(2^1  +  <*5)  »  ^  =  ®l  • 

The  mverse  of  the  fourth>order  tensor  (3.28)  is  (Kunin,  1983) 
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(3.30) 


6 

Me  =  • 

tsl 


The  scalar  parameters  Ci  are  related  to  parameters  6,,  defined  in  (3.29),  by 


I  -h.  -h.  1. 

A’fcs’66^’ 


where  A  =  i?  —  fcj  —  6|  +  . 

Returning  to  the  basis  T,  the  elastic  compliance  (3.30)  becomes 


€ 

At.  =  E  I' . 

tsl 


(3.31) 


(3.32) 


where  the  coefficients  d,-  are: 


dj  =  c«  ,  <**  =  -  CJ  -  Co) 

ds  =  2^”“®^  +  cj  —  2c3  —  204  +  cs)  ,  d#  =  +  cj  —  C3  +  C4  +  cg)  (3.33) 

^  *“  ^^®*  *”  ®*^  *  ds  —  +  C2  +  3c3  +  C4  —  4c5  +  ce) . 


In  the  eiq;tlidt  component  form,  (3.32)  reads 

Alyfci  =  di^iSikSji  +  SiiSjk)  +  d2  Stj^ki 

+  d3(%nfcni  +  nfnjikt)  % 

j  (3.34) 

+  d5-(tf,fcn,n,  +  ftifikSji  +  Siiitj-nk  +  iiiniSjk) 

+  ds  nitijUicni . 

In  dew  of  (3.29)  and  (3.31),  C4  -  C3  s  -2(c3  +  C4)  ;  hence,  from  (3.33),  dz  —  .  Thus, 

Me  satisfies  the  redprodty  relation  Miju  =  Muiij  ,  as  wdl  as  the  symmetry  properties 
Miiki  =  Mjikt  —  Mijik  ,  imposed  by  the  symmetry  of  the  stress  and  strain  tensors.  The  tensor 
(3.32),  or  (3.34),  is  the  compliance  tensor  of  the  considered  transversdy  isotropic  materiaL  The 
relationships  between  the  parameters  d,*  and  the  elastic  properties  of  the  transversely  isotropic 
material  are  easy  to  establish,  and  anil  be  derived  later  in  the  report  (Section  4.3). 
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3.5.  APPLICATIONS 


3.5.1.  Uniaxial  Loading 

Consider  first  a  prismatic  specimen  subjected  to  uniaxial  (tensile  or  compressive)  stress  a  direct¬ 
ed  along  the  axis  of  rotational  symmetry  (Fig.  3.6).  In  this  case  Oij  =  cr  n,nj  ,  where  n,  =  Sis  ■ 
The  corresponding  strain  is  Eij  =  [(1  +  t/)ninj  ~vSij]e  ,  where  v  (=  vis  =  vss)  is  the  Poisson’s 
coefficient  in  the  (1,2)  plane  of  elastic  symmetry,  while  a  and  e  are  the  longitudinal  stress  and 
strain.  Application  of  (3.22)  provides  in  this  case  the  following  stress-strain  relationship 

Oij  =  {A  n,nj  +  B  Sij)€  .  (3.35) 


r 


ia 

(a)  (b) 


Figure  S.6.  (a)  Planar  distribution  of  cracks:  all  cracks  have  their  normals  parallel  to  the  longi¬ 
tudinal  direction;  (b)  Cylindrical  distribution  of  cracks:  all  cracks  have  their  norm^  orthogonal 
to  the  longitudinal  direction. 


The  parameters 


i4  =  (1  -  Pn)[2^  (1  -  Pn)  -  ^0  Ap]  +  2|/(1  -  Pm)[P0  (1  “  Pm)  +  Aq  Ap  ]  (3.36) 


and 


B  =  (1  -  Pm)[Ao  (1  -  Pn)  -  2«/(Ao  +  /io)(l  -  pm)]  (3.37) 

serve  as  the  repositories  of  the  accumulated  damage.  Equating  (3.35)  with  <7,^-  =  o  riiiij  ,  it 
follows  that  B  =  0  and  A€  =  o’  .  The  condition  .F  =  0  provides  the  expression  for  the 
Poisson’s  coefficient 


1/  =  «/o  ,  (3.38) 

1  ~  Pm 

where  i/io  =  Ao/2(Ao  +  po)  i<  the  Poisson’s  coefficient  of  the  undamaged  isotropic  material. 

E^qpression  (3.38)  requires  further  explanation.  If  pm  =  0  ,  (3.38)  gives  i/  =  (1  —  pn)vo  , 
which  is  the  value  of  the  Poisson’s  coefficient  of  the  transversely  isotropic  material  containing 
cracks  embedded  in  parallel  planes  normal  to  the  longitudinal  axis  (Fig.  3  6a),  and  loaded  in 
tension.  The  compression  does  not  activate  the  cracks  in  Fig.  3.6a,  and  material  behaves 
in  compression  as  though  it  was  undamaged  (u  —  vio). 

If  s  0  ,  (3.38)  gives  uss(l  —  pm)~^*^  ,  which  is  the  Poisson’s  coefficient  of  the  transversely 
isotropic  material  with  a  qrlindiical  distribution  of  cracks  (Fig.  3.6b),  loaded  in  compression. 
This  crack  distribution  is  activated  only  in  the  presence  of  a  compressive  lonptudinal  stress, 
while  it  remaiiu  inactive  when  subjected  to  a  longitudinal  tensile  stress. 

PVom  the  remaining  condition,  iie  =  a  ,  it  itdlows  that  the  longitudinal  Young’s  modulus  is 
o/e  =  .4  .  In  view  of  (3.38)  and  identity  Aq  s  2i^(po  +  Aq)  ,  (3.36)  reduces  to 


il  =  (l-Pn)*^,  (3.39) 

where  JEb  =  2po(I  +  ^))  i>  the  Young’s  modulus  of  the  undamaged  isotropic  matrix.  When 
=  0  ,  the  lon^tudinal  Young’s  modulus  becomes  A  =  Bo  regardless  of  the  existence  of 
lonptudinal  cracks,  i.e.,  whether  p„  is  equal  to  aero  or  not.  H  p^  =  0  ,  the  volumetric  strain 
during  the  compression  test  is 


Bki 


1 


2^0 
~  Pm 


(3.40) 
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or,  by  introducing  the  initial  bulk  modulus  kq  3=  l?o/3(l  -  2i/o)  , 


Ekk  = 


21^0 

1  -  2vo 


Pm 

1  pm 


)a  . 


(3.41) 


During  a  gradual  increase  of  compression  a  (above  some  threshold  value),  Pm  increases  corre¬ 
spondingly.  In  the  case  when  an  appropriate  damage  law  =  Pm(<7)  is  additionally  provided 
(from  experimental  evidence  or  by  micromechanical  considerations),  (3.40)  completely  specifies 
the  volumetric  strain-stress  response.  A  qualitative  stress-strain  dependence,  depicted  in  Fig. 
3.7a,  replicates  the  basic  features  of  experimentally  observed  behavior  during  an  unconfined 
compression  test  (  tr  <  0  )  (see,  for  example,  Jaeger  and  Cook,  1976).  Volumetric  strain  is 
reduced  to  zero  at  deformation  states  for  which  the  crack  density  is  Pm  =  1  -  2vo  • 


o  o 


(a) 


(b) 


Figure  S.7.  (a)  Volumetric  stress-strain  response  corresponding  to  eqn  (3.41)  of  text:  c  is  the 
compressive  stress,  and  Ekk  i*  the  corresponding  volumetric  strain,  kq  is  the  undamaged  bulk 
modulus;  (b)  Stress-strain  response  corresponding  to  eqn  (3.49)  of  text:  a  is  the  magnitude  of 
the  biaxial  tension,  and  e  is  the  corresponding  strain,  Fb  aiid  ^  the  undamaged  Young’s 
modulus  and  Poisson’s  ratio. 
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3.5.2.  Biaxial  Loading 


Coniider  next  the  biaxial  tensile  loading  0\i  =  021  =  it  in  the  (1,2)  plane  of  elastic  sjrmmetry. 
The  corresponding  stress  and  str^  tensors  are: 

fTij  —  (iij  ~  Tiinj  )<T 

,  /  X  ,  (3-42) 

-  (1  +  • 

where  e  is  the  magnitude  of  strain  in  the  (1,2)  plane,  while  i>  =  2i^i/(l  -  i/u)  is  the  corre¬ 
sponding  (effective)  Poisson’s  coefficient.  From  (3.22)  it  further  follows  that 

Oij  =  (C  n.nj  +  D  f,j)«  ,  (3.43) 

in  which  the  accumulated  damage  is  recorded  using  the  parameters 

C  =  -2(1  -  Pm)[/io  (1  -  ^m)  +  Ao  Ap]  -  i>(l  -  p«)[2/io  (1  -  Pn)  -  Ao  Ap  ]  (3.44) 

and 


^  =  (1  “  Pm)l2(Ao  +  Po)  (1  -  Pm)  -  i>Ao  (1  -  Pn)] .  (3.45) 

Equating  the  first  of  equations  (3.42)  inth  (3.43)  it  foUows  that  P  =  -C  ,  and  therefore 

where  Pq  =  2Ao/(Ao+2/io)  —  2mi/(1— i^)  is  the  effective  Poisson’s  coefficient  of  the  undamaged 
material.  Substitution  of  (3.46)  into  (3.45)  then  pves 

D  =  (3«) 

i.  — 

The  stress-strain  equation  (3.43),  therefore,  becomes 

ffii  =  (1  -  Pm)’  i^ij  -  n,«,)e  .  (3.48) 

The  (o’,  e)  relation  can  now  be  deri'ved  from  (3.48)  in  conjunction  with  (3.42),  as 

»  =  (1  -  .  (3.«) 

The  (o,e)  response  is  completdy  specified  by  (3.49)  and  the  appropriate  damage  law  pm  = 
Pm(<r) .  A  qualitative  stress-strain  dependence  is  depicted  in  Fig.  3.7b. 
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4.0  MICROMECHANICALLY  INSPIRED  DAMAGE  MODELS 
FOR  BRITTLE  DEFORMATION  OF  ROCKS 

4.1.  ELASTIC  BODY  CONTAINING  A  PENNY-SHAPED  CRACK 

Expressions  for  the  elastic  stiffness  and  compliance  can  also  be  derived  on  the  basis  of  microme¬ 
chanical  models.  Consider  a  single  penny-shaped  crack  embedded  in  an  infinite  isotropic  elastic 
solid,  uniformly  loaded  at  infinity.  We  decompose  this  problem,  as  usual,  into  two  problems: 
that  of  the  body  without  a  crack,  loaded  at  infinity  and  that  of  the  body  with  a  crack 
appropriately  loaded  over  the  crack  faces  (*).  Correspondingly,  the  local  strain  e  can  be  written 
as  the  sum  of  the  strains  belonging  to  two  problems,  i.e.  e  =  e°  -I-  e*  .  The  averaged  strains 
(e  =  ■^  fyC  dV  as  V  — »  oo  )  are  decomposed  in  the  same  manner,  e  =  -{-  e*  .  Let  a-  be 

the  remote  loading  and  the  elastic  compliance  of  the  virgin  material  without  the  crack. 
Introducing  Aie  as  the  average  elastic  compliance  of  the  body  with  a  crack,  and  denoting  by 
Af  *  the  compliance  mapping  <r  to  e*  ,  e*  =  AC*  :  cr  ,  it  follows  that 

Me=Mj*  +  M*  .  (4.1) 

The  compliance  AC*  ,  being  the  Hessian  of  the  complementary  strain  energy  '9*  2^*  * 


8^9* 

r*"  be  conveniently  determined  by  observing  that  is  equal  to  the  energy  release  associated 
with  the  self-similar  crack  growth  from  aero  to  the  current  sixe  a.  As  shown  by  Budiansky  and 
Bice  (1973),  this  energy  can  be  expressed  as 


•'0  ® 

where  M  is  the  M-conservation  integral  of  fracture  mechanics.  The  M  integral  can  be  written 
in  terms  of  the  J  integral  by  means  of  a  line  integral  along  the  crack  edge  I 


M 


aJdl. 


Substituting  (4.4)  into  (4.3),  it  follows  that 


(4.4) 


(«) 
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In  the  dote  neighborhood  of  the  crack  edge,  the  stress  and  strain  states  are  a  combination  of 
plane  strain  and  antiplane  shear.  Thus,  the  energy  release  rate  or  the  J  integral  can  be  expressed 
in  terms  of  the  corresponding  stress-intensity  factors  Kj  (7  =  /,  II,  III)  as 

J  =  (^/  +  ^h)  +  2^  ^///  ■  (4.6) 

The  expression  (4.6)  can  be  conveniently  rewritten  as  (Sumarac,  1987) 


J  =  CijKiKj  , 


(4.7) 


where 


ClJ  =  [(1  -  Vo)6jj  -f-  ]  . 

Consequently,  substituting  (4.7)  into  (4.5)  gives 

*•=/“(/  CtjKjKj  dl)da.  (4.8) 

Jo  Ji 


Figure  4.1.  Penny-shaped  crack  of  radius  a  and  drcumference  1:  (l',2',3')  denotes  local  crack 
coordinate  system,  direction  1'  bang  coinddent  with  normal  to  crack  plane  m;  angles  9  and  ^ 
define  orientation  of  vector  m  relative  to  global  coodinate  system  (1,2,3). 
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For  ft  penny-shftped  crack  the  streis-inteniity  factors,  written  in  the  symmetricized  form,  are 
(lUft,  et  a/.,  1985): 


Kij  =  ((^12  +  <^2i)co8a  +  (<r{3  +  tr^i)BinQ  ]  (4.9) 

Kni  =  f  +  <^M)«ma  -  (trja  +  fr^,)co8Q  ]  , 

where  erlj  are  the  stress  components  in  the  crack  coordinate  system,  and  a  is  an  angle  defined 
in  Fig.  4.1.  The  stress  component  is  assumed  to  be  tensile.  Differentiating  (4.9)  it  follows 
that: 


=  ^(ira)'^^  — ^ —  [(^«i^i2  +  SiiSji)ca*a  +  {S^Sj^  +  fi3^ii)BU2a]  (4.10) 

m  • 


^<3  » 


~  —(to)'/*  ^ [(^n^i2  +^,-2tfyi)iina  -  {SnSfl  -f  dj3jj])cosa] . 


d(Tij  T 


2  —  uo 


If  0|]  is  a  compzes^e  streH,  JT/  s  0  along  with  the  right>haad  side  of  the  first  expression  in 
(4.10). 

The  components  of  the  compliance  tensor  M*  in  the  local  (crack)  coordinate  system  are 


....  '  »  r  ilr^  BKjdKj  . 

**>>»  -  d<r!j»ci,  “  *4  «crt,  ^  ' 


(4.11) 


Since 


BKjdKj  _l-vo  .BKidKi  dKndKn  1  BKujdKui 

2/10  ^Bc•^jBa•^^  Bc\^  B<ry^  2ito  Be^ij  Bc’^i 


(4.12) 


substitution  of  (4.10)  and  (4.12)  into  (4.11)  leads  to  the  fcdlowing  expression  for  the  compliance 
tensor 
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(4.13) 


[(^tl^iJ  +  Si2i}l)ifkl^l2  +  SkzSn)  +  (SiiSjs  +  tft3^jl)(^kl^/3  +  ^fc3^n)]}  • 

The  Heaviside  step  function  H{eT[i)  is  introduced  in  (4.13)  to  simultaneously  treat  the  possi¬ 
bility  of  tensile  and  compressive  stress  components  a'l^.  In  the  local  crack  coordinate  system, 
the  normal  to  the  crack  plane  has  the  components 

mj  =  5.1  .  (4.14) 

The  expression  for  the  compliance  tensor  (4.13)  thus  becomes 

Mi/ki  - 

(4.15) 

The  expression  (4.15)  can  be  farther  rearranged  as 


M*iju  =  y 

j  (4.16) 

The  required  symmetry  properties  M^jh  =  Mj/u  =  M*/i^  —  Mhfj  clearly  hold.  Using  the 
fourth-order  tensors  forming  the  integrity  basis  1%  introduced  in  Section  3.3,  expression  (4.16) 
can  be  rewritten  in  compact  form  as 


161-»b  1 


^  “  *^)^(‘^«)  -  2]  iy/J?'  }  •  (4.17) 

I«  («7).  ifjS'  and  are  the  tensors  defined  by  (3.17),  in  terms  of  the  components  of 

the  normal  m  relative  to  the  local  crack  coordinate  system  (  m(-  =  Su  ).  If  cri,  £  0  .  (4.17) 
further  reduces  to 


Z  2- VO  2n^ 
used  by  Krajdnovic  and  Fandla  (1986). 


,  .  _  16  1  -  »ID  1  3  /O  rS  m'  _  _  76  m'\ 


(4.18) 
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4.2.  AVERAGING  PROCEDURE 


The  compliance  tensor  expressed  in  (4.17)  relative  to  the  local  crack  coordinate  system,  can  be 
written  in  the  global  coordinate  system  using  the  coordinate  transformation 


~  QiaQjfi  QhyQlS  > 

where  Q  is  the  orthogonal  tensor  of  the  transformation  between  the  two  coordinate  systems. 
Substitution  of  (4.17)  gives 


{  2  I?iS  + 1(2  -  n,)j(m’'  ■  .r  •  m)  -  2]  IfjS  }  ,  (4.19) 

where: 


~  +  mimkSji  +  Siimjmk  +  m,m/fjt) 

^jki  =  rmmjmkmi 


(4.20) 


are  the  fonrth>order  tensors  combining  the  Kronecker  delta  and  the  unit  vector  m,  which  is 
expressed  in  the  global  coordinate  system  as 


m  3s  {cos^  COS0  >cos^  sin9  }  • 

If  ‘  O’  <  m  >  0  (o-  is  the  stress  tensor  with  components  in  the  global  coordinate 

system),  (4.19)  reduces  to 

Expressions  (4.21)  and  (4.18)  are  equal,  except  that  in  (4.18)  the  components  of  the  normal  m 
are  expressed  in  the  local  crack  coordinate  system,  while  in  (4.21)  they  are  expressed  relative  to 
the  global  comdinate  system.  The  superscript  m  indicates  that  the  reference  is  to  the  plane  of 
the  crack  (Fig.  4.1). 

Consider  now  the  case  of  many  cracks.  U  all  cracks  have  the  same  normal  m,  the  average 
compliance  (neglecting  the  direct  crack  interaction)  is  JA*  =  NM.*  ,  where  N  is  the  number 
of  cracks  per  unit  volume.  Next,  let  the  crack  distribution  be  such  that  aU  normals  to  the  crack 
planes  have  the  same  an^  ^  —  constant.  N^lecting  direct  interaction  between  adjacent  cracks 
(dilute  distribution  uf  cracks),  the  average  compliance  is 
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dB  = 


(4.22) 


M 


•  -nr 

“  2»  Jq 

jy  16 1  -  1 

2t  3  2  -  Wo  2/<o 


dS  -vq 


dB  )  . 


In  (4.22)  the  itress  state  is  assumed  to  be  such  that  •  o’  •  m  >  0  for  all  m,  allowing 
substitution  of  (4.21)  for  the  components  of  the  compliance  tensor  . 

In  general,  for  an  arbitrary  angle  ^ 


m,m,  =  SixSjxtxM^^CM^B  +  f.jfjjcos’ Rein’d  + 

+  +  drtdji)cos’^iintfcos^  +  (BiiSjs  +  Si3Sji)tm^os^osff  (4.23) 

+  (fi2^j3  +  di3^i2)«ia^«^«a^  • 


Consider  first  an  important  special  case  for  which  ^  =  0.  The  normal  m  has  components 
{cosd  ,  sind  ,  0}  .  In  this  case  (4.23)  reduces  to 


miiTij  =  fiidji  cos^d  +  SaSji  ain*d  +  (dadj-s  +  SaSji)  sindcosd  .  (4.24) 


Hence 


»a» 

/  mimj  dB  =  ir(dAdji  +  dijdia)  -  »(d«i  -  dodp) 

Jfk 


(4-25) 


Introduce  the  vector  n,  normal  to  m,  having  the  components 


=  do 


(4.26) 


in  the  global  coordinate  qrstem  .  The  integral  (4.25)  can  be  accordin^y  written  as 


inch  that 


mimj  dd  =  w(dy  -  tiiUj) , 


(4.27) 


J  l!iSM=  +  lutik) 

—  (dffcn^ni  +  ninitSji  +  SunjUk  +  njntd^-fc)] . 


(4.28) 


XJnng  the  fourth-order  tensors  of  the  int^rity  basis  I*  introduced  in  (3.17),  the  integr<.I  (4*28) 
can  be  rewritten  in  compact  form  as 


40 


f2it 

/  I*"*  =  T(r  -  I*")  . 

Jo 

(4.29) 

To  ev&luate  the  integral 

.2*  ,2ir 

/  Ifjll}  d0=  rmmjmkmi  d0  , 

•'0  •'0 

(4.30) 

expressions  of  the  type  (4.24)  are  substituted  for  products  m,m^  and  m^mf.  Performing  the 
intergration,  one  obtains 

1  ^jkl  +  ^•2^i2^fc2^/2) 

•'0  ’ 

+  (f«l^j2^i:l^/2  +  ^il^i2^i2^n  ^f2^il^*3^/2 

+  +  ^»l^il^k2^/2  +  ^t2^i2^kl^/l)]  • 

(4.31) 

After  a  somewhat  delicate  rearrangement  of  terms,  (4.31)  can  be  cast  in  the  following  remarkably 
simple  form 

^  d$  s  ^[(fy  -  n.n,)(^w  -  n*n,) 

+  i^jk  —  njnk){Sii  —  «»nj)  +  (^w  -  nkni)(Sji  —  njiti)] , 

(4.32) 

or 

IT 

/  -^jw  +  ^ti^Jk)  +  —  {Sijnifni  +  tiiiijSki) 

-  (fiknjtii  +  +  n,n,f^fc)  +  3(nin,n*nj)]  . 

(4.33) 

Utiliang  the  fourth-order  tensors  forming  the  integrity  basis  (3.17),  (4.33)  can  be  written  as 

I*"*  <W  =  ?^(2  I* -H’ - 1®**  - -  4  I®** -1- 3  I®*)  . 

Jo  ’ 

(4.34) 

finally,  the  average  compliance  due  to  the  considered  crack  distribution  is  derived  by  substituting 
(4.28)  and  (4.33)  into  (4.22) 

•^yw  =  ^  (^ik^Ji  +  ^ii^jk)  -  Vo  +  VQ  (^yn/fenj  -I-  ninjSki) 

-  (2  -  »b)  (^*fc»i«i  +  ninkSji  +  -b  iwnjd,*)  -  3i^  (n,n,n^nj)] . 

(4.35) 
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The  nondimensional  scalar  quantity  b>  =  Na^,  originally  introduced  by  Budiansky  and  O’Connell 
(1976),  represents  a  micromechanical  damage  parameter  (measure)  defining  the  density  of  the 
considered  crack  distribution  within  the  representative  volume  element.  Note  that  M’  is  linearly 
proportional  to  ut,  i.e.  /A*  =  uM*  ,  where  M*  is  the  constant  tensor  given  by 


.  21-1^  1 


M’  =  r 


[  2(4  -  i>b)  -  I/O  +  I/O  (I^"  +  I^")  -  4(2  -  i/o)  I*"  -  3vo  ] 


3  2-1/0  2/io 

Finally,  if  the  crack  distribution  is  dilute  and  isotropic,  the  average  compliance  is 


'0  •'-ir/2 

By  using  (4.20)  and  (4.23)  it  can  bo  shown  that: 


M*  =  —  /  A4*  cos^  d(f> 

4ir  Jq  •/_»/2 


de 


(4.36) 


/  /  I*"  cos^  ^  I* 

Jo  J.r/7  3 

fir  fr/2  4 

/  /  !«"*  COS^  ^  (2  +  I*)  . 

•/O  -f-w/i  15 

The  average  compliance  is  derived  by  substituting  (4.37)  into  (4.21)  and  (4.36) 


(4J7) 


It  is  easily  shown  that  (4.38)  produces  the  compliance  components  identical  to  those  fpven  by 
Equation  (20)  of  Hoxii  and  Nemat-Nasser  (1983). 


4.3.  EFFECTIVE  COMPLIANCE  TENSOR 

The  cmnpliance  tensor  of  the  undamaged,  isotropic  and  homogeneous  elastic  matrix  is 

M°iih  =  ^  tiifk,  ] . 

Using  the  basis  I*  introduced  in  Section  3.3,  above  can  be  rewritten  as 


A4.®  =  ^(  V  ) 

2/io'  l  +  m 


(4.39) 


The  overall  (effective)  compliance  is  derived  by  superposing  (4.39)  and  the  expression  for  the 
average  compliance  JA*.  U  the  crack  distribution  is  isotropic,  i.e.  if  JA*  is  defined  by  (4.38), 
the  overall  compliance  becomes 
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1/ 


(4.40) 


=  ;^(  I‘  -  7-^^  )  , 

2/i'‘  1  +  1/ 


2/i'  1  + 

where  the  damage-dependent  shear  modulus  and  Poisson’s  ratio  are  given  by: 


M  = 


V  = 


45(2- I/O  ) 


fio 


45(2  —  i/o)  +  32(1  —  i/o)(5  —  vo)u> 
45(2  -  i/to)  +  16(1  -  i/^)w 
45(2  -  i/o)  +  16(1  -  i/^)(10  -  3i/o)a; 


(4.41) 


*^0  • 


If  the  crack  distribution  is  of  the  type  shown  in  Fig.  3.6b,  such  that  (4.35)  applies,  it  follows 
that 


-  (2  -  Wo)  {Siknjni  +  ninkSjt  +  Sunjiik  +  n,n,tf,jt)  -  3i/o  (n<n_,n*nj)]} 


(4.42) 


Expression  (4.42)  represents  the  compliance  tensor  of  a  transversely  isotropic  material  whose 
plane  of  symmetry  is  normal  to  the  direction  n.  Indeed,  (4.42)  can  be  rewritten  as 

6 

Ai.  =  2  Ci  r  ,  (4.43) 

«sl 

with  the  obvioiu  eoqpressions  for  the  parameters  Ci.  In  terms  of  the  elastic  moduli  and  Poisson’s 
ratios  in  two  orthogonal  directions,  the  parameters  Ci  are: 

V 

~E 

Vi/ 

E>  (4.44) 

1  +  2^  w  .  1  1 

E'  ~  E'^2ji~  IP  ' 

Here,  E  is  the  Young’s  modulus  in  the  plane  of  isotropy  and  E'  is  in  the  direction  normal  to 
it.  Also,  V  is  the  Poisson’s  ratio  characterizing  transverse  contraction  in  the  plane  of  isotropy 
when  tension  is  applied  in  the  same  plane,  whOe  i/  is  the  Poisson’s  ratio  obtained  when  tension 
is  applied  normal  to  the  plane  of  isotropy.  is  the  shear  modulus  for  any  plane  perpendicular 
to  the  plane  of  isotropy  (Lekhnitskii,  1981).  These  five  material  parameters  can  be  derived  from 
(4.44)  in  terms  of  the  parameters  C,-  as: 


Ci  = 


l  +  v 


1 

2/s  ’ 


(75  =  1-1, 

II*  II 


Cj 

Cs 

C, 
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(4.45) 


E  = 


f?'  = 


C,  +  Cj  ’ 


Cl  +  Cj  +  2C3  +  Cj  +  Cg  ’ 


i/  = 


V  = 


m'  = 


2C1  +  C5 


C2 

Cl  +  C2  ’ 

Cj  +  C3 

Cl  4-  C2  +  2C3  +  C5  +  Ce 


4.4.  EFFECTIVE  STIFFNESS  TENSOR 

To  derive  the  expression  for  the  stiffness  tensor,  being  the  inverse  of  the  compliance  tensor  (4.42) 
from  Section  4.3,  it  is  convenient  to  nse  the  fourth-order  tensors  of  the  basis  J'  ,  defined  by  (3.27). 
In  this  basis,  the  tensor  (4.35)  takes  the  form 

M'  =  0, 1  (2  -  «,)  (J>"  -  J»")  +  2  J‘"  +  (4  -  ►»)  J»"  I .  (4.46) 

Since: 


I»  =  i(3J»-J»-3J®  +  J^), 

the  initial  compliance  ^veu  by  (4.39),  can  be  expressed  in  terms  of  the  J*  basis.  The 
overall  compliance  M,  s  M,^  +  M*  ,  consequently  becomes 


A4.= 


(4.47) 


where: 


_  2  ”  Mo  ,  ^  \ 

‘■=2(l  +  «,)  +  3<‘”"“>" 
—  3*4) 

^  2(1  + Mo)  * 

,  .  8(1 -Mb) 

'‘  =  ‘+3(2T^)"- 


C2  = 


*4) 


C4  =  - 


2(1  +  Mo)  3 
*4) 


-  7(1  -  *4))  w  . 


2(1  + *4))  ’ 


(4.48) 


,^.4<l-«j)(4^__ 
3(2 -Mo) 


The  tensor  representation  (4.47)  has  the  explicit  inverse 


£.  =  A4.”'  =  2aio  Y;biJ\ 

ml 


(4.49) 
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where 


{h, . = 


(4.50) 


and  A  =  cj  -  c|  -  c5  +  cj  .  Using  the  transformation  rule  (3.27),  (4.49)  can  be  expressed  in 
terms  of  the  original  I*  basis  as 


S  tti  I’ 


The  coefficients  Oi  are  related  to  coefficients  6^  by: 


(4.51) 


®1  =  ^6  1 


=  2(^1  -  ^2  -  h) 


as  =  ^(— hi  +  hi  —  2h3  —  2h4  +  he)  ,  04  =  x(~hi  +  hs  —  ha  +  h4  +  he)  (4.52) 

2  « 


Os  =  2(hs  —  hfi)  , 


06  =  ^(3hi  +  ha  +  3h3  +  h4  —  4hs  +  he) 

A 


which  is  identical  to  the  relationship  (3.33)  existing  between  the  coefficients  d,  and  Cj.  In  the 
component  form,  (4.51)  reads 


=  2/io  [  +  ai  htihfc/ 

+  os|(h<fcnyn/  +  niiikSji  +  h.in^n*  +  n,n/h,k) 
+  06  ninjtifcni  ]  . 


(4.53) 


Eaqirestion  (4.53)  is  the  exact  inverse  of  (4.42).  If  the  damage  parameter  u  is  snffidently  small 
so  th»-t  the  quadratic  and  higher-order  terms  in  uf  can  be  n^lected,  coefficients  of  (4.53)  simplify 
accordingjly  and  become: 


as  =  I  (1  -  I'o)  w 


(4.54) 


45 


06  =  - 


2 

3 


«^o(l  -  t^) 
(2-i^)(l-2^,) 


(13  —  2vo)  u)  . 


It  is  easy  to  prove  that  for  this  case  the  components  of  (4.53)  axe  identical  to  the  components  of 
the  stillness  tensor  obtwed  by  Nemat-Nasser  and  Hori  (1990)  (their  Equation  (5.11b)).  Also, 
from  (4.53)  and  (4.54)  it  clearly  follows  that  Ce  =  +  (>;  L*  ,  where 


L*  =  [  (oi  -  1)  +  (02  —  ^  ^  + 1^)  +  Os  I®  +  06  I®  ]  ,  (4.55) 

and  r*®  =  2/io(r  +  I*)  . 


4.5.  AXIAL  COMPRESSION  OF  A  LATERALLY  CONFINED  SPECIMEN 

Consider  a  cylindrical  specimen  subjected  to  an  axial  compression  a  +  p  and  lateral  confining 
pressure  p.  The  stress  difference  (between  axial  smd  lateral  stress),  commonly  used  in  rock 
mechanics,  is  a.  The  state  of  macroscopic  stress  is  therefore 


=  “(P  ^0  +  o  n,n,) ,  (4.56) 

where  n,-  s  is  a  unit  vector  coUineax  with  the  longitudinal  axis  of  the  cylinder.  Even 
though  the  applied  tractions  are  compressive,  local  tensile  stresses  may  arise  in  the  vicinity  of  the 
microstructural  inhomogeneities  (pores,  pre-existing  cracks,  rigid  inclusions,  etc.).  Introducing 
the  influence  tensor  B  of  the  given  microstructural  inhomogendty  (Hiil,  1967),  one  can  write 
the  local  stress  tensor  as  <r*  =  B  :  tr  ^  where  B  depends  on  the  location  and  topology  of  the 
particular  defect.  At  a  low  level  of  lateral  confinement  the  local  tensile  stresses  may  suffice  to 
nucleate  new  microcracks  and  propagate  easting  ones.  The  ensuing  deformation  is  brittle,  i.e., 
the  inelastic  deformation  is  directly  related  to  the  formation  of  new  internal  surfaces  (nudeation 
of  new  cracks  and  growth  of  ones)  in  the  specimen.  Assuming  that  these  microcracks 

develop  in  planes  with  normals  perpendicular  to  the  longitudinal  (axial)  directicm  n  (Fig.  3.6b), 
the  initially  isotropic  material  becomes  transversdy  isotropic.  In  the  simplest  representation, 
the  local  tensile  stress  driving  the  microcrack  growth  (Fig.  4.2),  can  be  written  as 

"ti  =  (<«  -  >w>i)  .  («  w) 

where  the  magnitude  of  the  local  tension  o*  depends  on  the  microstructure,  i.e.,  on  the  defects 
such  as  porosity  or  void  distributions,  and  the  stress  difference  a.  The  influence  tensor  B 
corresponding  to  the  local  stress  pven  by  the  simplified  representation  (4.57)  is  not  unique. 
One  of  its  possible  representatioiu,  satisfying  <r*  =  F  :  o’ ,  is  B  =  a[I^  —  ^(1  —  I^)]  , 

where  ot  =  a*  fa  and  —  aff  . 
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The  straia  reaulting  from  the  crack  opening!  due  to  the  streu  (4.57)  is 


(4.58) 

where  is  the  compliance  tensor,  given  by  (4.35).  Hence, 

«*i  =  |(1  -  Wo)  w  ^  {Sij  -  nitij)  .  (4.59) 

The  strain  due  to  macroscopic  stress  in  the  material  without  damage  (i.e.,  with  the  initial,  frozen 
microstructure)  is 


,  (4.60) 

i.e.,  in  view  of  (4.56)  and  (4.39), 

Defining  the  total  strain  as  €ij  s  i  it  follows  that  the  longitudinal  strain  is 

e33  =  ~(<r  +  (l-2^,)pl,  (4.62) 

where  Eo  =  2^(1  +  Mi)  u  ^he  initial  Young’s  modulus  of  the  undamaged  material.  The  total 
volumetric  strain  is 


hi  (4.63),  kq  a:  £b/3(l  -  2mo)  is  the  bulk  modulus  cS  the  undamaged  matrix.  Consequently, 
as  a  result  of  the  introduced  amplifications,  the  longitudinal  strain  (4.62)  k  not  influenced  by 
the  damage.  This  is  consistent  with  the  assumption  that  the  components  of  the  crack-induced 
displacement  discontinuities  in  the  Erection  parallel  to  the  specimen  axis  are  zero.  On  the  other 
hand,  the  non-zero  strain  components  =  <22  contribute  to  damage-affected  volumetric 
strain,  pven  by  (4.63).  With  the  additionally  provided  relationship  defining  the  variation  of 
the  parameter  ua*  as  a  function  of  increasing  stress  a  and  pressure  p,  (4.63)  suffices  for  the 
determination  of  the  volumetric  strain-stress  response.  The  qualitative  dependence  is  as  shown 
in  Fig.  3.7a. 
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If  the  ipecimen  if  not  confined  (p  =  0),  the  Poitton’i  coefficient  i/'  in  the  (1,2)  plane  of  elastic 
symmetry  due  to  compression  in  the  direction  of  the  X3  axis  (  1/'  =  -(11/(33  =  -(22/^33  = 
;(1  -  fkk/tss)  )  follows  from  (4.62)  and  (4.63) 


t/'  =  Mo  +  |(1  -  vi)ut—  .  (4.64) 

Again,  if  the  relationship  between  and  the  applied  stress  is  known,  (4.64)  specifies  the  change 

of  the  Poisson’s  ratio  t/'  caused  by  continuing  d^radation,  defined  by  the  damage  parameter  u>. 

I” 


Figure  4.2.  Local  tennle  streu  a*  (which  driyet  the  microcrack  growth)  arises  in  the  vicinity 
of  the  pore. 
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5.0  SUMMARY  AND  CONCLUSIONS 


The  objective  of  this  report  was  to  itunmarixe  on-going  research  directed  towards  the  formulation 
of  a  rational  continuum  model  of  a  porous  rock  containing  a  large  number  of  microcracks.  The 
ultimate  goal  is  to  provide  a  versatile  physically-based  model  suitable  for  large-scale,  finite 
element-type  computations.  There  is  a  considerable  body  of  evidence  that  a  variety  of  important 
features  of  rock  deformation  (such  as  dilatancy,  compaction,  brittle-ductile  transition,  etc.) 
cannot  be  modeled  without  consideration  of  the  rock  fabric  and  its  influence  on  micromechanical 
processes.  On  the  other  hand,  purely  micromechanical  models  are  computationally  inefficient 
and,  therefore,  not  suitable  for  such  purposes. 

In  view  of  the  program  objective,  the  only  rational  solution  is  to  develop  a  micromechanically- 
inspired  constitutive  model.  The  proposed  modd  will  retain  the  form  and  framework  of  con¬ 
ventional  phenomenolopcal  theories.  The  details  of  the  model,  however,  such  as  the  selection 
of  mathematical  representations  of  the  damage  parameters,  forms  of  evolution  laws,  etc.,  will 
be  constructed  in  concert  with  micromechanical  considerations  of  micromechanisms  observed 
experimentally  in  real  porous  rocks  (sandstones  and  limestones).  Consequently,  the  material 
parameters  of  the  model  will  be  experimentally  identifiable.  However,  their  numerical  values 
will  require  some  minimal  amount  of  curve  fitting.  An  illustration  of  micromechanical  modeling 
was  presented  in  Sections  2.2  and  2.3. 

The  present  study  is  at  the  stage  where  it  focuses  on  brittle  response  by  examining  the  ge¬ 
ometry  of  microcracks  and  their  influence  on  the  response.  Sections  3.1  and  3.2  concentrated 
on  the  physical  description  of  damage  defined  by  frictional  crack  surfaces  in  planes  of  various 
orientations  passing  through  a  material  p<unt.  Various  approximations  of  a  typical  damage 
distribution  in  porous  rocks  were  illustrated  in  Sections  3.1  and  3.2.  It  is  unlikely  in  real  ap¬ 
plications  that  a  detailed  crack  density  diitiibution  will  ever  be  available.  Thus,  it  was  argued 
that  the  second-order  tensor  damage  parameter  approximates  the  expected  crack  density  distri¬ 
bution  with  sufficient  accuracy  in  the  case  of  proportional  loading.  The  case  of  non-proportional 
loading  is  scheduled  to  be  addressed  in  the  sequel  to  this  study. 

Rigorous,  novel,  and  elegant  derivation  of  the  effective  (overall  or  equivalent)  stiffness  and  ccnn- 
pliance  tensors  was  presented  in  Sections  3.3  and  3.4.  The  formulas  for  the  exact  inversions 
of  these  tensors  were  derived  to  enhance  intended  computational  efficiency.  Particular  cases  of 
isotropic  and  transversdy  isotropic  damage  distributions  were  derived  not  only  as  an  illustration 
but  more  importantly  as  a  guide  to  the  determination  of  material  parameters  and  the  damage 
evolution  law.  Applications  to  the  ciue  of  uniaxial  and  biajdal  stress  (compressive  and/or  tensile) 
loadings  were  illustrated  in  Section  3.5. 

Section  4.0  provided  additional  guidelines  in  rdating  the  proposed  model  to  the  physics  of  the 
process  on  the  microscale.  The  emphasis  was  again  on  the  enhancement  of  tensor  manipulations 
of  the  sti&esses  and  compliances.  At  each  step  the  model  constants  were  rdated  to  measurable 
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m&terial  parameters  (see  expressions  (4.45),  for  example).  The  example  in  Section  4.5  provided 
an  illustration  of  how  the  proposed  model  can  be  used  in  estimating  accumulated  damage  during 
a  given  loading  program.  In  addition  to  acoustic  emission  tests,  this  will  provide  a  valuable  means 
for  determining  damage  evolution  laws  and/or  the  kinematics  of  the  “damage  surface”. 

In  summary,  the  present  model  provides  a  general  framework  for  a  rational,  micromechanically- 
inspired,  continuum  model  for  a  porous  rock.  The  model  provides  a  novel  viewpoint  of  the  entire 
class  of  damage  models.  For  the  first  time,  clear  arguments  have  been  provided  to  facilitate 
the  selection  of  the  simplest  (yet  sufficiently  accurate)  representation  of  damage.  Microcrack 
distributions,  damage  meastires  (parameters),  and  material  parameters  are  firmly  related  to 
each  other  allowing  for  unambiguous  experimental  identifications  of  the  quantities  needed  for 
modeling. 

At  this  point  the  analysis  is  limited  to  purdy  brittle  deformation  processes  and  proportion¬ 
al  loadings.  Incorporation  of  ductile  effects  and  the  transition  to  rate  models  (necessary  for 
non-proportional  loadings)  will  be  addressed  in  the  sequel  to  this  study  and  when  requisite  ex¬ 
perimental  data  become  available.  In  order  to  define  the  damage  evolution  law(s)  and,  perhaps, 
define  a  “damage  surface”,  it  will  be  necessary  to  carefully  design  an  experimental  program 
combining  measurements  of  elastic  moduli  in  three  orthogonal  directions  along  with  acoustic 
emissions  measurements.  These  experiments  will  have  to  be  performed  at  several  confining 
pressure  levels  for  different  loading  programs  in  order  to  develop  a  reliable  basis  for  the  deter¬ 
mination  of  the  evolution  law(s).  A  major  problem,  of  course,  is  our  inability  to  measure  the 
microcrack  density  at  each  point  directly.  Thus,  it  becomes  necessary  to  relate  the  microcrack 
distribution  to  macrostiffiiesses,  as  was  done  in  this  report,  and  infer  the  former  from  the  latter. 
The  importance  of  this  task  is  not  evident  from  the  problems  considered  in  this  report.  However, 
in  the  case  of  non-proportional  loading  it  will  be  necessary  to  know  the  orientation  of  cracks  in 
order  to  estimate  the  discontinuous  changes  in  elastic  moduli  associated  with  the  changes  in  the 
signs  of  the  normal  stresses. 

A  few  simple  examples  were  used  to  illustrate  the  remarkable  versatility  and  efficiency  of  the 
model  in  replicating  the  salient  trends  of  the  conudered  phenomena.  Even  in  its  present  state  of 
development  the  proposed  model  allows  for  a  umple  identification  of  the  material  parameters  and 
their  experimental  measurement.  In  particular,  it  appears  possible  to  determine  the  approximate 
distribution  of  damage  densities  by  measuring  the  components  of  the  stiffness  tensor.  The 
indenter  test  proposed  by  Zarka  and  IVdat  (1977)  could  be  particularly  appropriate  for  this 
task. 

The  most  important  conclusion  derived  from  the  work  reported  heron  is  that  it  is  indeed  pos¬ 
sible  to  complete  this  task  and  formulate  a  constitutive  theory  for  rocks  which  will  satisfy  all 
requirements  of  efficiency  and  accuracy.  The  general  form  of  this  theory  has  already  been  put 
together  in  this  report.  However,  a  significant  effort  must  stiU  be  mounted  in  order  to  achieve 
the  ambitious  goal  of  this  overall  research  project. 
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